# =============================================================================
# FILE 4. MAIN ANALYSES
# =============================================================================
# PURPOSE:
#   Runs the primary statistical analyses, including
#   data loading/merging and a suite of regression models and diagnostics.
#
# INPUTS (expected in working directory):
#   - cf_flanders_2023_data.rds
#   - cf_flanders_2024_data.rds
#   - cf_flanders_data.rds
#   - pf_flanders_2023_data.rds
#   - pf_flanders_2024_data.rds
#   - pf_flanders_data.rds
# =============================================================================

# =============================================================================
# 1. WORKING DIRECTORY + FILE LOCATIONS
# =============================================================================

getwd()

# Set project directory (update it to your own directory that contains the input files)
setwd("updated_path_here")

# =============================================================================
# 2. Create Function to get modal values
# =============================================================================

getmode <- function(v) {
  uniqv <- unique(v)
  uniqv[which.max(tabulate(match(v, uniqv)))]
}

# =============================================================================
# 3. Create party labels
# =============================================================================

party.labels <- c("Other", "Preferred Party")
party.labels.pvda <- c("Other", "PVDA")
party.labels.groen <- c("Other", "Groen")
party.labels.vooruit <- c("Other", "Vooruit")
party.labels.ovld <- c("Other", "Open Vld")
party.labels.cdv <- c("Other", "CD&V")
party.labels.nva <- c("Other", "N-VA")
party.labels.vb <- c("Other", "Vlaams Belang")

# =============================================================================
# 4. Load data
# =============================================================================

# Load citizen data
# --- Data import: 2023
cf.flanders.2023.data <- readRDS("cf_flanders_2023_data.rds")
# --- Data import: 2024
cf.flanders.2024.data <- readRDS("cf_flanders_2024_data.rds")
# --- Data import: 2023-2024
cf.flanders.data <- readRDS("cf_flanders_data.rds")

# Load politician data
# --- Data import: 2023
pf.flanders.2023.data <- readRDS("pf_flanders_2023_data.rds")
# --- Data import: 2024
pf.flanders.2024.data <- readRDS("pf_flanders_2024_data.rds")
# --- Data import: 2023-2024
pf.flanders.data <- readRDS("pf_flanders_data.rds")

# =============================================================================
# 5. MERGE citizen and politician data
# =============================================================================

# Rename vote column to party in citizen dataframe
cf.flanders.data <- cf.flanders.data %>%
  dplyr::rename(party = vote)

# Rename binary vote intention columns in citizen dataframe
cf.flanders.data <- cf.flanders.data %>%
  dplyr::rename(party_pvda = vote_pvda,
                party_groen = vote_groen,
                party_vooruit = vote_vooruit,
                party_ovld = vote_ovld,
                party_cdv = vote_cdv,
                party_nva = vote_nva,
                party_vb = vote_vb)

# Rename binary party affiliation columns in politician dataframe
pf.flanders.data <- pf.flanders.data %>%
  dplyr::rename(party_pvda = affiliation_pvda,
                party_groen = affiliation_groen,
                party_vooruit = affiliation_vooruit,
                party_ovld = affiliation_ovld,
                party_cdv = affiliation_cdv,
                party_nva = affiliation_nva,
                party_vb = affiliation_vb)

# Merge citizen and politician dataframes
flanders.data <- dplyr::bind_rows(cf.flanders.data, pf.flanders.data)

# Define variable type
flanders.data$party <- as.character(flanders.data$party)
flanders.data$party_pvda <- as.factor(flanders.data$party_pvda)
flanders.data$party_groen <- as.factor(flanders.data$party_groen)
flanders.data$party_vooruit <- as.factor(flanders.data$party_vooruit)
flanders.data$party_ovld <- as.factor(flanders.data$party_ovld)
flanders.data$party_cdv <- as.factor(flanders.data$party_cdv)
flanders.data$party_nva <- as.factor(flanders.data$party_nva)
flanders.data$party_vb <- as.factor(flanders.data$party_vb)
flanders.data$interest <- as.factor(flanders.data$interest)
flanders.data$sex <- as.factor(flanders.data$sex)
flanders.data$age <- as.numeric(flanders.data$age)
flanders.data$ideology <- as.numeric(flanders.data$ideology)
flanders.data$survey <- as.factor(flanders.data$survey)
flanders.data$type <- as.factor(flanders.data$type)
flanders.data$seats_pvda <- as.factor(flanders.data$seats_pvda)
flanders.data$seats_groen <- as.factor(flanders.data$seats_groen)
flanders.data$seats_vooruit <- as.factor(flanders.data$seats_vooruit)
flanders.data$seats_cdv <- as.factor(flanders.data$seats_cdv)
flanders.data$seats_ovld <- as.factor(flanders.data$seats_ovld)
flanders.data$seats_nva <- as.factor(flanders.data$seats_nva)
flanders.data$seats_vb <- as.factor(flanders.data$seats_vb)
flanders.data$govt_pvda <- as.numeric(flanders.data$govt_pvda)
flanders.data$govt_groen <- as.numeric(flanders.data$govt_groen)
flanders.data$govt_vooruit <- as.numeric(flanders.data$govt_vooruit)
flanders.data$govt_cdv <- as.numeric(flanders.data$govt_cdv)
flanders.data$govt_ovld <- as.numeric(flanders.data$govt_ovld)
flanders.data$govt_nva <- as.numeric(flanders.data$govt_nva)
flanders.data$govt_vb <- as.numeric(flanders.data$govt_vb)
flanders.data$seats_accuracy <- as.numeric(flanders.data$seats_accuracy)
flanders.data$sophistication_index <- as.numeric(flanders.data$sophistication_index)

# Rescale cumulative Brier score from 0 to 1
flanders.data <- flanders.data %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  mutate(govt_brier = (govt_brier - min(govt_brier, na.rm = TRUE)) / (max(govt_brier, na.rm = TRUE) - min(govt_brier, na.rm = TRUE)))

# Save dataframe in Stata format
write_dta(flanders.data, "flanders_data.dta")

# Save dataframe in R format
saveRDS(flanders.data, "flanders_data.rds")

# =============================================================================
# 6. Percentage of correct seat forecasts
# =============================================================================

# Results reported in section 4.2 of the paper (seat accuracy)

seat_party_map <- tibble::tibble(
  seats_col   = c("seats_pvda_correct",
                  "seats_groen_correct",
                  "seats_vooruit_correct",
                  "seats_cdv_correct",
                  "seats_ovld_correct",
                  "seats_nva_correct",
                  "seats_vb_correct"),
  party_match = c("PVDA",
                  "Groen",
                  "Vooruit",
                  "CD&V",
                  "Open Vld",
                  "N-VA",
                  "Vlaams Belang")
)

# Base long data (one row per respondent x seats_col)
base_long <- flanders.data %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  dplyr::select(survey, type, party, dplyr::all_of(seat_party_map$seats_col)) %>%
  tidyr::pivot_longer(
    cols = dplyr::all_of(seat_party_map$seats_col),
    names_to = "seats_col",
    values_to = "correct"
  ) %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  dplyr::left_join(seat_party_map, by = "seats_col") %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  dplyr::mutate(
    correct = as.numeric(correct),
    type_label = dplyr::if_else(type == 1, "politicians", "citizens")
  ) %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  dplyr::filter(!is.na(correct))  # always drop NA correctness values

summarise_version <- function(df) {
  results_long <- df %>%
    # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
    dplyr::mutate(
      match_group = dplyr::if_else(!is.na(party) & party == party_match,
                                   "party == match",
                                   "party != match")
    ) %>%
    dplyr::group_by(survey, type_label, seats_col, party_match, match_group) %>%
    # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
    dplyr::summarise(
      n = dplyr::n(),
      n_correct = sum(correct == 1),
      pct_correct = 100 * mean(correct == 1),
      .groups = "drop"
    )
  
  results_wide <- results_long %>%
    # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
    dplyr::select(survey, type_label, seats_col, party_match, match_group, n, n_correct, pct_correct) %>%
    tidyr::pivot_wider(
      names_from  = match_group,
      values_from = c(n, n_correct, pct_correct)
    ) %>%
    dplyr::arrange(survey, type_label, seats_col)
  
  list(long = results_long, wide = results_wide)
}

# A) Treat party NA as non-match (included)
res_including_party_na <- summarise_version(base_long)

# B) Filter out party NA
# --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
res_excluding_party_na <- summarise_version(dplyr::filter(base_long, !is.na(party)))

out_path <- "flanders_seats_accuracy_by_party_type_survey.xlsx"
writexl::write_xlsx(
  list(
    "wide_including_partyNA" = res_including_party_na$wide,
    "long_including_partyNA" = res_including_party_na$long,
    "wide_excluding_partyNA" = res_excluding_party_na$wide,
    "long_excluding_partyNA" = res_excluding_party_na$long
  ),
  path = out_path
)

out_path

# =============================================================================
# 7. Regression models: Expectations
# =============================================================================

# Fit ordered logistic regression models with triple interaction: seat expectations
model.seats_pvda <- clm(seats_pvda ~ as.factor(party_pvda)*as.factor(type)*as.factor(survey) + sex + age, data = flanders.data)
model.seats_groen <- clm(seats_groen ~ as.factor(party_groen)*as.factor(type)*as.factor(survey) + sex + age, data = flanders.data)
model.seats_vooruit <- clm(seats_vooruit ~ as.factor(party_vooruit)*as.factor(type)*as.factor(survey) + sex + age, data = flanders.data)
model.seats_cdv <- clm(seats_cdv ~ as.factor(party_cdv)*as.factor(type)*as.factor(survey) + sex + age, data = flanders.data)
model.seats_ovld <- clm(seats_ovld ~ as.factor(party_ovld)*as.factor(type)*as.factor(survey) + sex + age, data = flanders.data)
model.seats_nva <- clm(seats_nva ~ as.factor(party_nva)*as.factor(type)*as.factor(survey) + sex + age, data = flanders.data)
model.seats_vb <- clm(seats_vb ~ as.factor(party_vb)*as.factor(type)*as.factor(survey) + sex + age, data = flanders.data)

# Fit OLS regression models with triple interaction: government expectations
model.govt_pvda <- lm(govt_pvda ~ as.factor(party_pvda)*as.factor(type)*as.factor(survey) + sex + age, data = flanders.data)
model.govt_groen <- lm(govt_groen ~ as.factor(party_groen)*as.factor(type)*as.factor(survey) + sex + age, data = flanders.data)
model.govt_vooruit <- lm(govt_vooruit ~ as.factor(party_vooruit)*as.factor(type)*as.factor(survey) + sex + age, data = flanders.data)
model.govt_cdv <- lm(govt_cdv ~ as.factor(party_cdv)*as.factor(type)*as.factor(survey) + sex + age, data = flanders.data)
model.govt_ovld <- lm(govt_ovld ~ as.factor(party_ovld)*as.factor(type)*as.factor(survey) + sex + age, data = flanders.data)
model.govt_nva <- lm(govt_nva ~ as.factor(party_nva)*as.factor(type)*as.factor(survey) + sex + age, data = flanders.data)
model.govt_vb <- lm(govt_vb ~ as.factor(party_vb)*as.factor(type)*as.factor(survey) + sex + age, data = flanders.data)

# =============================================================================
# 8. Regression model: Correct seat forecasts (PVDA)
# =============================================================================

# Convert the dependent variable to a factor
flanders.data$seats_pvda_correct <- as.factor(flanders.data$seats_pvda_correct)

# Fit a cumulative link model (ordinal logistic regression)
model.seats_pvda_correct <- clm(seats_pvda_correct ~ as.factor(party_pvda)*as.factor(type)*as.factor(survey) + 
                                  sex + age, data = flanders.data)

# Display model summary without the correlation matrix between coefficients
summary(model.seats_pvda_correct, correlation = FALSE)

# Compute estimated marginal means (EMMs) of 'type' within each level of 'party_pvda'
# This shows the first difference: differences in predicted values across 'type'
pvda.int1 <- emmeans(model.seats_pvda_correct, ~ type | party_pvda * survey)

# Compute the second difference (i.e., difference in first differences)
pvda.int2 <- pairs(pairs(emmeans::regrid(pvda.int1)), by = NULL)
pvda.int2  # View the second differences

# =============================================================================
# 9. Regression model: Correct seat forecasts (Groen)
# =============================================================================

# Convert the dependent variable to a factor
flanders.data$seats_groen_correct <- as.factor(flanders.data$seats_groen_correct)

# Fit a cumulative link model (ordinal logistic regression)
model.seats_groen_correct <- clm(seats_groen_correct ~ as.factor(party_groen)*as.factor(type)*as.factor(survey) + 
                                   sex + age, data = flanders.data)

# Display model summary without the correlation matrix between coefficients
summary(model.seats_groen_correct, correlation = FALSE)

# Compute estimated marginal means (EMMs) of 'type' within each level of 'party_groen'
# This shows the first difference: differences in predicted values across 'type'
groen.int1 <- emmeans(model.seats_groen_correct, ~ type | party_groen * survey)

# Compute the second difference (i.e., difference in first differences)
groen.int2 <- pairs(pairs(emmeans::regrid(groen.int1)), by = NULL)
groen.int2  # View the second differences

# =============================================================================
# 10. Regression model: Correct seat forecasts (Vooruit)
# =============================================================================

# Convert the dependent variable to a factor
flanders.data$seats_vooruit_correct <- as.factor(flanders.data$seats_vooruit_correct)

# Fit a cumulative link model (ordinal logistic regression)
model.seats_vooruit_correct <- clm(seats_vooruit_correct ~ as.factor(party_vooruit)*as.factor(type)*as.factor(survey) + 
                                     sex + age, data = flanders.data)

# Display model summary without the correlation matrix between coefficients
summary(model.seats_vooruit_correct, correlation = FALSE)

# Compute estimated marginal means (EMMs) of 'type' within each level of 'party_vooruit'
# This shows the first difference: differences in predicted values across 'type'
vooruit.int1 <- emmeans(model.seats_vooruit_correct, ~ type | party_vooruit * survey)

# Compute the second difference (i.e., difference in first differences)
vooruit.int2 <- pairs(pairs(emmeans::regrid(vooruit.int1)), by = NULL)
vooruit.int2  # View the second differences

# =============================================================================
# 11. Regression model: Correct seat forecasts (CD&V)
# =============================================================================

# Convert the dependent variable to a factor
flanders.data$seats_cdv_correct <- as.factor(flanders.data$seats_cdv_correct)

# Fit a cumulative link model (ordinal logistic regression)
model.seats_cdv_correct <- clm(seats_cdv_correct ~ as.factor(party_cdv)*as.factor(type)*as.factor(survey) + 
                                 sex + age, data = flanders.data)

# Display model summary without the correlation matrix between coefficients
summary(model.seats_cdv_correct, correlation = FALSE)

# Compute estimated marginal means (EMMs) of 'type' within each level of 'party_cdv'
# This shows the first difference: differences in predicted values across 'type'
cdv.int1 <- emmeans(model.seats_cdv_correct, ~ type | party_cdv * survey)

# Compute the second difference (i.e., difference in first differences)
cdv.int2 <- pairs(pairs(emmeans::regrid(cdv.int1)), by = NULL)
cdv.int2  # View the second differences

# =============================================================================
# 12. Regression model: Correct seat forecasts (Open Vld)
# =============================================================================

# Convert the dependent variable to a factor
flanders.data$seats_ovld_correct <- as.factor(flanders.data$seats_ovld_correct)

# Fit a cumulative link model (ordinal logistic regression)
model.seats_ovld_correct <- clm(seats_ovld_correct ~ as.factor(party_ovld)*as.factor(type)*as.factor(survey) + 
                                  sex + age, data = flanders.data)

# Display model summary without the correlation matrix between coefficients
summary(model.seats_ovld_correct, correlation = FALSE)

# Compute estimated marginal means (EMMs) of 'type' within each level of 'party_ovld'
# This shows the first difference: differences in predicted values across 'type'
ovld.int1 <- emmeans(model.seats_ovld_correct, ~ type | party_ovld * survey)

# Compute the second difference (i.e., difference in first differences)
ovld.int2 <- pairs(pairs(emmeans::regrid(ovld.int1)), by = NULL)
ovld.int2  # View the second differences

# =============================================================================
# 13. Regression model: Correct seat forecasts (N-VA)
# =============================================================================

# Convert the dependent variable to a factor
flanders.data$seats_nva_correct <- as.factor(flanders.data$seats_nva_correct)

# Fit a cumulative link model (ordinal logistic regression)
model.seats_nva_correct <- clm(seats_nva_correct ~ as.factor(party_nva)*as.factor(type)*as.factor(survey) + 
                                 sex + age, data = flanders.data)

# Display model summary without the correlation matrix between coefficients
summary(model.seats_nva_correct, correlation = FALSE)

# Compute estimated marginal means (EMMs) of 'type' within each level of 'party_nva'
# This shows the first difference: differences in predicted values across 'type'
nva.int1 <- emmeans(model.seats_nva_correct, ~ type | party_nva * survey)

# Compute the second difference (i.e., difference in first differences)
nva.int2 <- pairs(pairs(emmeans::regrid(nva.int1)), by = NULL)
nva.int2  # View the second differences

# =============================================================================
# 14. Regression model: Correct seat forecasts (Vlaams Belang)
# =============================================================================

# Convert the dependent variable to a factor
flanders.data$seats_vb_correct <- as.factor(flanders.data$seats_vb_correct)

# Fit a cumulative link model (ordinal logistic regression)
model.seats_vb_correct <- clm(seats_vb_correct ~ as.factor(party_vb)*as.factor(type)*as.factor(survey) + 
                                sex + age, data = flanders.data)

# Display model summary without the correlation matrix between coefficients
summary(model.seats_vb_correct, correlation = FALSE)

# Compute estimated marginal means (EMMs) of 'type' within each level of 'party_vb'
# This shows the first difference: differences in predicted values across 'type'
vb.int1 <- emmeans(model.seats_vb_correct, ~ type | party_vb * survey)

# Compute the second difference (i.e., difference in first differences)
vb.int2 <- pairs(pairs(emmeans::regrid(vb.int1)), by = NULL)
vb.int2  # View the second differences

# =============================================================================
# 15. Interaction: Correct seat forecasts (PVDA)
# =============================================================================

# Fit logistic regression model
model.seats_pvda_correct <- glm(seats_pvda_correct ~ party_pvda*type*survey + 
                                  sex + age, family="binomial", data = flanders.data)

# Display model summary
summary(model.seats_pvda_correct)

# Predicted probabilities of correct forecast by type and party preference
# Not computed: Nearly all PVDA respondents predicted more seats for PVDA

# Predicted probabilities of correct forecast by party preference only
pvda.pred.party.prop <- ggeffect(model.seats_pvda_correct, terms = c("party_pvda"))
pvda.pred.party.ref  <- ggpredict(model.seats_pvda_correct, terms = c("party_pvda"))

# Predicted probabilities of correct forecast by respondent type only (citizen or politician)
pvda.pred.type.prop <- ggeffect(model.seats_pvda_correct, terms = c("type"))
pvda.pred.type.ref  <- ggpredict(model.seats_pvda_correct, terms = c("type"))

# Compute Average Marginal Effects (AMEs) for the interaction between party preference and type
pvda.eff.int <- effect(term="party_pvda*type*survey", mod=model.seats_pvda_correct)

# Convert the effect object to a data frame for easier manipulation and plotting
pvda.eff.int.data <- as.data.frame(pvda.eff.int)

# Add a column labeling the party for clarity in plots/tables
pvda.eff.int.data$party <- "PVDA (+5)"

# Show the effect data frame (predicted probabilities and confidence intervals)
pvda.eff.int.data

# Test significance of the interaction between party preference and type using estimated marginal means (EMMs)
pvda.int.f <- emmeans(model.seats_pvda_correct, ~ type | party_pvda * survey)
pvda.int.s <- pairs(pairs(emmeans::regrid(pvda.int.f)), by = NULL)

# Plot the first difference effects with comparison indicators
plot(pvda.int.f, comparisons = TRUE)

# Prepare factor levels for plotting to ensure correct order on x-axis and legend
pvda.eff.int.data$type <- factor(pvda.eff.int.data$type, levels = c("0", "1"))
pvda.eff.int.data$party_pvda <- factor(pvda.eff.int.data$party_pvda, levels = c("0", "1"))

# Save the interaction plot to a high-resolution PNG file
png("vb_correct_int_g.png", units="in", width=6, height=4, res=1200)

# Create a ggplot showing predicted probabilities of correct forecast by type and party preference
# --- Plotting: figure creation starts here. Review aesthetics (aes) and any facet/scales for readability.
pvda.eff.int.g <- ggplot(pvda.eff.int.data, aes(x=type, y=fit, color=party_pvda, group=party_pvda)) +
  geom_line(size=.5, show.legend = TRUE, aes(linetype=party_pvda, color=party_pvda)) +      # Lines for each party group
  geom_point(size=.5, show.legend = FALSE) +                                             # Points for predicted values
  geom_errorbar(aes(ymin=fit-se, ymax=fit+se, color=party_pvda), width = 0.05, show.legend = FALSE) +  # Confidence intervals
  labs(title = "PVDA") +                                                                # Plot title
  scale_x_discrete(labels = c("0" = "Citizens", "1" = "Politicians")) +                  # Label x-axis categories
  scale_y_continuous(name = "Predicted Probability", limits=c(0, 1),                      # y-axis scale and label
                     breaks = seq(0, 1, .2)) +
  facet_wrap(~ survey, ncol = 2) +
  theme_minimal() +                                                                      # Minimal theme for clean look
  theme(axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9),
        axis.title = element_text(size = 9), 
        axis.title.y = element_text(margin = margin(t = 0, r = 5, b = 0, l = 0)),
        axis.title.x = element_blank(),
        legend.title = element_text(size = 8), legend.text=element_text(size=8),
        legend.key.size = unit(.8, 'cm'),
        plot.title = element_blank(),
        plot.margin=unit(c(0,.01,.3,0), "cm"), legend.position="bottom") +
  scale_colour_manual(values=c("#636363", "#990000"), name="Party Preference",            # Custom colors for party groups
                      labels=c("0"="Other","1"="PVDA"), guide = guide_legend(reverse=TRUE)) +
  scale_fill_manual(values=c("#636363", "#990000"), name="Party Preference",              # Fill colors (if needed)
                    labels=c("0"="Other","1"="PVDA"), guide = guide_legend(reverse=TRUE)) +
  scale_linetype_manual(values=c("longdash","solid"), name="Party Preference",            # Line types for party groups
                        labels=c("0"="Other","1"="PVDA"), guide = guide_legend(reverse=TRUE))

# Print the plot object to the file
pvda.eff.int.g

# Close the PNG device to save the file
dev.off()

# =============================================================================
# 16. Interaction: Correct seat forecasts (Groen)
# =============================================================================

# Fit logistic regression model
model.seats_groen_correct <- glm(seats_groen_correct ~ party_groen*type*survey + 
                                   sex + age, family="binomial", data = flanders.data)

# Display model summary
summary(model.seats_groen_correct)

# Predicted probabilities of correct forecast by type and party preference
groen.pred.int.prop <- ggeffect(model.seats_groen_correct, terms = c("type", "party_groen", "survey"))
groen.pred.int.ref <- ggpredict(model.seats_groen_correct, terms = c("type", "party_groen", "survey"))

# Predicted probabilities of correct forecast by party preference only
groen.pred.party.prop <- ggeffect(model.seats_groen_correct, terms = c("party_groen"))
groen.pred.party.ref <- ggpredict(model.seats_groen_correct, terms = c("party_groen"))

# Predicted probabilities of correct forecast by respondent type only (citizen or politician)
groen.pred.type.prop <- ggeffect(model.seats_groen_correct, terms = c("type"))
groen.pred.type.ref <- ggpredict(model.seats_groen_correct, terms = c("type"))

# Compute Average Marginal Effects (AMEs) for the interaction between party preference and type
groen.eff.int <- effect(term="party_groen*type*survey", mod=model.seats_groen_correct)

# Convert the effect object to a data frame for easier manipulation and plotting
groen.eff.int.data <- as.data.frame(groen.eff.int)

# Add a column labeling the party for clarity in plots/tables
groen.eff.int.data$party <- "Groen (–5)"

# Show the effect data frame (predicted probabilities and confidence intervals)
groen.eff.int.data

# Test significance of the interaction between party preference and type using estimated marginal means (EMMs)
groen.int.f <- emmeans(model.seats_groen_correct, ~ type | party_groen * survey)
groen.int.s <- pairs(pairs(emmeans::regrid(groen.int.f)), by = NULL)

# Plot the first difference effects with comparison indicators
plot(groen.int.f, comparisons = TRUE)

# Prepare factor levels for plotting to ensure correct order on x-axis and legend
groen.eff.int.data$type <- factor(groen.eff.int.data$type, levels = c("0", "1"))
groen.eff.int.data$party_groen <- factor(groen.eff.int.data$party_groen, levels = c("0", "1"))

# Save the interaction plot to a high-resolution PNG file
png("groen_correct_int_g.png", units="in", width=6, height=4, res=1200)

# Create a ggplot showing predicted probabilities of correct forecast by type and party preference
# --- Plotting: figure creation starts here. Review aesthetics (aes) and any facet/scales for readability.
groen.eff.int.g <- ggplot(groen.eff.int.data, aes(x=type, y=fit, color=party_groen, group=party_groen)) +
  geom_line(size=.5, show.legend = TRUE, aes(linetype=party_groen, color=party_groen)) +
  geom_point(size=.5, show.legend = FALSE) +
  geom_errorbar(aes(ymin=fit-se, ymax=fit+se, color=party_groen), width = 0.05, show.legend = FALSE) +
  labs(title = "Groen") +
  scale_x_discrete(labels = c("0" = "Citizens", "1" = "Politicians")) +
  scale_y_continuous(name = "Predicted Probability", limits=c(0, 1),
                     breaks = seq(0, 1, .2)) +
  facet_wrap(~ survey, ncol = 2) +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9),
        axis.title = element_text(size = 9), 
        axis.title.y = element_text(margin = margin(t = 0, r = 5, b = 0, l = 0)),
        axis.title.x = element_blank(),
        legend.title = element_text(size = 8), legend.text=element_text(size=8),
        legend.key.size = unit(.8, 'cm'),
        plot.title = element_blank(),
        plot.margin=unit(c(0,.01,.3,0), "cm"), legend.position="bottom") +
  scale_colour_manual(values=c("#636363", "#04847A"), name="Party Preference",
                      labels=c("0"="Other","1"="Groen"), guide = guide_legend(reverse=TRUE)) +
  scale_fill_manual(values=c("#636363", "#04847A"), name="Party Preference",
                    labels=c("0"="Other","1"="Groen"), guide = guide_legend(reverse=TRUE)) +
  scale_linetype_manual(values=c("longdash","solid"), name="Party Preference",
                        labels=c("0"="Other","1"="Groen"), guide = guide_legend(reverse=TRUE))

# Print the plot object to the file
groen.eff.int.g

# Close the PNG device to save the file
dev.off()

# =============================================================================
# 17. Interaction: Correct seat forecasts (Vooruit)
# =============================================================================

# Fit logistic regression model
model.seats_vooruit_correct <- glm(seats_vooruit_correct ~ party_vooruit*type*survey + 
                                     sex + age, family="binomial", data = flanders.data)

# Display model summary
summary(model.seats_vooruit_correct)

# Predicted probabilities of correct forecast by type and party preference
vooruit.pred.int.prop <- ggeffect(model.seats_vooruit_correct, terms = c("type", "party_vooruit", "survey"))
vooruit.pred.int.ref <- ggpredict(model.seats_vooruit_correct, terms = c("type", "party_vooruit", "survey"))

# Predicted probabilities of correct forecast by party preference only
vooruit.pred.party.prop <- ggeffect(model.seats_vooruit_correct, terms = c("party_vooruit"))
vooruit.pred.party.ref <- ggpredict(model.seats_vooruit_correct, terms = c("party_vooruit"))

# Predicted probabilities of correct forecast by respondent type only (citizen or politician)
vooruit.pred.type.prop <- ggeffect(model.seats_vooruit_correct, terms = c("type"))
vooruit.pred.type.ref <- ggpredict(model.seats_vooruit_correct, terms = c("type"))

# Compute Average Marginal Effects (AMEs) for the interaction between party preference and type
vooruit.eff.int <- effect(term="party_vooruit*type*survey", mod=model.seats_vooruit_correct)

# Convert the effect object to a data frame for easier manipulation and plotting
vooruit.eff.int.data <- as.data.frame(vooruit.eff.int)

# Add a column labeling the party for clarity in plots/tables
vooruit.eff.int.data$party <- "Vooruit (+5)"

# Show the effect data frame (predicted probabilities and confidence intervals)
vooruit.eff.int.data

# Test significance of the interaction between party preference and type using estimated marginal means (EMMs)
vooruit.int.f <- emmeans(model.seats_vooruit_correct, ~ type | party_vooruit * survey)
vooruit.int.s <- pairs(pairs(emmeans::regrid(vooruit.int.f)), by = NULL)

# Plot the first difference effects with comparison indicators
plot(vooruit.int.f, comparisons = TRUE)

# Prepare factor levels for plotting to ensure correct order on x-axis and legend
vooruit.eff.int.data$type <- factor(vooruit.eff.int.data$type, levels = c("0", "1"))
vooruit.eff.int.data$party_vooruit <- factor(vooruit.eff.int.data$party_vooruit, levels = c("0", "1"))

# Save the interaction plot to a high-resolution PNG file
png("vooruit_correct_int_g.png", units="in", width=6, height=4, res=1200)

# Create a ggplot showing predicted probabilities of correct forecast by type and party preference
# --- Plotting: figure creation starts here. Review aesthetics (aes) and any facet/scales for readability.
vooruit.eff.int.g <- ggplot(vooruit.eff.int.data, aes(x=type, y=fit, color=party_vooruit, group=party_vooruit)) +
  geom_line(size=.5, show.legend = TRUE, aes(linetype=party_vooruit, color=party_vooruit)) +
  geom_point(size=.5, show.legend = FALSE) +
  geom_errorbar(aes(ymin=fit-se, ymax=fit+se, color=party_vooruit), width = 0.05, show.legend = FALSE) +
  labs(title = "Vooruit") +
  scale_x_discrete(labels = c("0" = "Citizens", "1" = "Politicians")) +
  scale_y_continuous(name = "Predicted Probability", limits=c(0, 1),
                     breaks = seq(0, 1, .2)) +
  facet_wrap(~ survey, ncol = 2) +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9),
        axis.title = element_text(size = 9), 
        axis.title.y = element_text(margin = margin(t = 0, r = 5, b = 0, l = 0)),
        axis.title.x = element_blank(),
        legend.title = element_text(size = 8), legend.text=element_text(size=8),
        legend.key.size = unit(.8, 'cm'),
        plot.title = element_blank(),
        plot.margin=unit(c(0,.01,.3,0), "cm"), legend.position="bottom") +
  scale_colour_manual(values=c("#636363", "#FF0000"), name="Party Preference",
                      labels=c("0"="Other","1"="Vooruit"), guide = guide_legend(reverse=TRUE)) +
  scale_fill_manual(values=c("#636363", "#FF0000"), name="Party Preference",
                    labels=c("0"="Other","1"="Vooruit"), guide = guide_legend(reverse=TRUE)) +
  scale_linetype_manual(values=c("longdash","solid"), name="Party Preference",
                        labels=c("0"="Other","1"="Vooruit"), guide = guide_legend(reverse=TRUE))

# Print the plot object to the file
vooruit.eff.int.g

# Close the PNG device to save the file
dev.off()

# =============================================================================
# 18. Interaction: Correct seat forecasts (CD&V)
# =============================================================================

# Fit logistic regression model
model.seats_cdv_correct <- glm(seats_cdv_correct ~ party_cdv*type*survey + 
                                 sex + age, family="binomial", data = flanders.data)

# Display model summary
summary(model.seats_cdv_correct)

# Predicted probabilities of correct forecast by type and party preference
cdv.pred.int.prop <- ggeffect(model.seats_cdv_correct, terms = c("type", "party_cdv", "survey"))
cdv.pred.int.ref <- ggpredict(model.seats_cdv_correct, terms = c("type", "party_cdv", "survey"))

# Predicted probabilities of correct forecast by party preference only
cdv.pred.party.prop <- ggeffect(model.seats_cdv_correct, terms = c("party_cdv"))
cdv.pred.party.ref <- ggpredict(model.seats_cdv_correct, terms = c("party_cdv"))

# Predicted probabilities of correct forecast by respondent type only (citizen or politician)
cdv.pred.type.prop <- ggeffect(model.seats_cdv_correct, terms = c("type"))
cdv.pred.type.ref <- ggpredict(model.seats_cdv_correct, terms = c("type"))

# Compute Average Marginal Effects (AMEs) for the interaction between party preference and type
cdv.eff.int <- effect(term="party_cdv*type*survey", mod=model.seats_cdv_correct)

# Convert the effect object to a data frame for easier manipulation and plotting
cdv.eff.int.data <- as.data.frame(cdv.eff.int)

# Add a column labeling the party for clarity in plots/tables
cdv.eff.int.data$party <- "CD&V (–3)"

# Show the effect data frame (predicted probabilities and confidence intervals)
cdv.eff.int.data

# Test significance of the interaction between party preference and type using estimated marginal means (EMMs)
cdv.int.f <- emmeans(model.seats_cdv_correct, ~ type | party_cdv * survey)
cdv.int.s <- pairs(pairs(emmeans::regrid(cdv.int.f)), by = NULL)

# Plot the first difference effects with comparison indicators
plot(cdv.int.f, comparisons = TRUE)

# Prepare factor levels for plotting to ensure correct order on x-axis and legend
cdv.eff.int.data$type <- factor(cdv.eff.int.data$type, levels = c("0", "1"))
cdv.eff.int.data$party_cdv <- factor(cdv.eff.int.data$party_cdv, levels = c("0", "1"))

# Save the interaction plot to a high-resolution PNG file
png("cdv_correct_int_g.png", units="in", width=6, height=4, res=1200)

# Create a ggplot showing predicted probabilities of correct forecast by type and party preference
# --- Plotting: figure creation starts here. Review aesthetics (aes) and any facet/scales for readability.
cdv.eff.int.g <- ggplot(cdv.eff.int.data, aes(x=type, y=fit, color=party_cdv, group=party_cdv)) +
  geom_line(size=.5, show.legend = TRUE, aes(linetype=party_cdv, color=party_cdv)) +
  geom_point(size=.5, show.legend = FALSE) +
  geom_errorbar(aes(ymin=fit-se, ymax=fit+se, color=party_cdv), width = 0.05, show.legend = FALSE) +
  labs(title = "CD&V") +
  scale_x_discrete(labels = c("0" = "Citizens", "1" = "Politicians")) +
  scale_y_continuous(name = "Predicted Probability", limits=c(0, 1),
                     breaks = seq(0, 1, .2)) +
  facet_wrap(~ survey, ncol = 2) +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9),
        axis.title = element_text(size = 9), 
        axis.title.y = element_text(margin = margin(t = 0, r = 5, b = 0, l = 0)),
        axis.title.x = element_blank(),
        legend.title = element_text(size = 8), legend.text=element_text(size=8),
        legend.key.size = unit(.8, 'cm'),
        plot.title = element_blank(),
        plot.margin=unit(c(0,.01,.3,0), "cm"), legend.position="bottom") +
  scale_colour_manual(values=c("#636363", "#FF6000"), name="Party Preference",
                      labels=c("0"="Other","1"="CD&V"), guide = guide_legend(reverse=TRUE)) +
  scale_fill_manual(values=c("#636363", "#FF6000"), name="Party Preference",
                    labels=c("0"="Other","1"="CD&V"), guide = guide_legend(reverse=TRUE)) +
  scale_linetype_manual(values=c("longdash","solid"), name="Party Preference",
                        labels=c("0"="Other","1"="CD&V"), guide = guide_legend(reverse=TRUE))

# Print the plot object to the file
cdv.eff.int.g

# Close the PNG device to save the file
dev.off()

# =============================================================================
# 19. Interaction: Correct seat forecasts (OPEN VLD)
# =============================================================================

# Fit logistic regression model
model.seats_ovld_correct <- glm(seats_ovld_correct ~ party_ovld*type*survey + 
                                  sex + age, family="binomial", data = flanders.data)

# Display model summary
summary(model.seats_ovld_correct)

# Predicted probabilities of correct forecast by type and party preference
ovld.pred.int.prop <- ggeffect(model.seats_ovld_correct, terms = c("type", "party_ovld", "survey"))
ovld.pred.int.ref <- ggpredict(model.seats_ovld_correct, terms = c("type", "party_ovld", "survey"))

# Predicted probabilities of correct forecast by party preference only
ovld.pred.party.prop <- ggeffect(model.seats_ovld_correct, terms = c("party_ovld"))
ovld.pred.party.ref <- ggpredict(model.seats_ovld_correct, terms = c("party_ovld"))

# Predicted probabilities of correct forecast by respondent type only (citizen or politician)
ovld.pred.type.prop <- ggeffect(model.seats_ovld_correct, terms = c("type"))
ovld.pred.type.ref <- ggpredict(model.seats_ovld_correct, terms = c("type"))

# Compute Average Marginal Effects (AMEs) for the interaction between party preference and type
ovld.eff.int <- effect(term="party_ovld*type*survey", mod=model.seats_ovld_correct)

# Convert the effect object to a data frame for easier manipulation and plotting
ovld.eff.int.data <- as.data.frame(ovld.eff.int)

# Add a column labeling the party for clarity in plots/tables
ovld.eff.int.data$party <- "Open Vld (–7)"

# Show the effect data frame (predicted probabilities and confidence intervals)
ovld.eff.int.data

# Test significance of the interaction between party preference and type using estimated marginal means (EMMs)
ovld.int.f <- emmeans(model.seats_ovld_correct, ~ type | party_ovld * survey)
ovld.int.s <- pairs(pairs(emmeans::regrid(ovld.int.f)), by = NULL)

# Plot the first difference effects with comparison indicators
plot(ovld.int.f, comparisons = TRUE)

# Prepare factor levels for plotting to ensure correct order on x-axis and legend
ovld.eff.int.data$type <- factor(ovld.eff.int.data$type, levels = c("0", "1"))
ovld.eff.int.data$party_ovld <- factor(ovld.eff.int.data$party_ovld, levels = c("0", "1"))

# Save the interaction plot to a high-resolution PNG file
png("ovld_correct_int_g.png", units="in", width=6, height=4, res=1200)

# Create a ggplot showing predicted probabilities of correct forecast by type and party preference
# --- Plotting: figure creation starts here. Review aesthetics (aes) and any facet/scales for readability.
ovld.eff.int.g <- ggplot(ovld.eff.int.data, aes(x=type, y=fit, color=party_ovld, group=party_ovld)) +
  geom_line(size=.5, show.legend = TRUE, aes(linetype=party_ovld, color=party_ovld)) +
  geom_point(size=.5, show.legend = FALSE) +
  geom_errorbar(aes(ymin=fit-se, ymax=fit+se, color=party_ovld), width = 0.05, show.legend = FALSE) +
  labs(title = "Open Vld") +
  scale_x_discrete(labels = c("0" = "Citizens", "1" = "Politicians")) +
  scale_y_continuous(name = "Predicted Probability", limits=c(0, 1),
                     breaks = seq(0, 1, .2)) +
  facet_wrap(~ survey, ncol = 2) +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9),
        axis.title = element_text(size = 9), 
        axis.title.y = element_text(margin = margin(t = 0, r = 5, b = 0, l = 0)),
        axis.title.x = element_blank(),
        legend.title = element_text(size = 8), legend.text=element_text(size=8),
        legend.key.size = unit(.8, 'cm'),
        plot.title = element_blank(),
        plot.margin=unit(c(0,.01,.3,0), "cm"), legend.position="bottom") +
  scale_colour_manual(values=c("#636363", "#005DAA"), name="Party Preference",
                      labels=c("0"="Other","1"="Open Vld"), guide = guide_legend(reverse=TRUE)) +
  scale_fill_manual(values=c("#636363", "#005DAA"), name="Party Preference",
                    labels=c("0"="Other","1"="Open Vld"), guide = guide_legend(reverse=TRUE)) +
  scale_linetype_manual(values=c("longdash","solid"), name="Party Preference",
                        labels=c("0"="Other","1"="Open Vld"), guide = guide_legend(reverse=TRUE))

# Print the plot object to the file
ovld.eff.int.g

# Close the PNG device to save the file
dev.off()

# =============================================================================
# 20. Interaction: Correct seat forecasts (N-VA)
# =============================================================================

# Fit logistic regression model
model.seats_nva_correct <- glm(seats_nva_correct ~ party_nva*type*survey + 
                                 sex + age, family="binomial", data = flanders.data)

# Display model summary
summary(model.seats_nva_correct)

# Predicted probabilities of correct forecast by type and party preference
# 'ggeffect' computes predicted marginal effects, setting other categorical variables to their proportions
nva.pred.int.prop <- ggeffect(model.seats_nva_correct, terms = c("type", "party_nva", "survey"))
# 'ggpredict' computes predicted probabilities with other categorical variables held at their reference levels
nva.pred.int.ref <- ggpredict(model.seats_nva_correct, terms = c("type", "party_nva", "survey"))

# Predicted probabilities of correct forecast by party preference only
nva.pred.party.prop <- ggeffect(model.seats_nva_correct, terms = c("party_nva"))
nva.pred.party.ref <- ggpredict(model.seats_nva_correct, terms = c("party_nva"))

# Predicted probabilities of correct forecast by respondent type only (citizen or politician)
nva.pred.type.prop <- ggeffect(model.seats_nva_correct, terms = c("type"))
nva.pred.type.ref <- ggpredict(model.seats_nva_correct, terms = c("type"))

# Compute Average Marginal Effects (AMEs) for the interaction between party preference and type
nva.eff.int <- effect(term="party_nva*type*survey", mod=model.seats_nva_correct)

# Convert the effect object to a data frame for easier manipulation and plotting
nva.eff.int.data <- as.data.frame(nva.eff.int)

# Add a column labeling the party for clarity in plots/tables
nva.eff.int.data$party <- "N-VA (–4)"

# Show the effect data frame (predicted probabilities and confidence intervals)
nva.eff.int.data

# Test significance of the interaction between party preference and type using estimated marginal means (EMMs)
# Calculate first differences: effect of 'type' within each level of 'party_nva'
nva.int.f <- emmeans(model.seats_nva_correct, ~ type | party_nva * survey)
# Calculate second differences: differences between first differences (interaction effects)
nva.int.s <- pairs(pairs(emmeans::regrid(nva.int.f)), by = NULL)

# Plot the first difference effects with comparison indicators
plot(nva.int.f, comparisons = TRUE)

# Prepare factor levels for plotting to ensure correct order on x-axis and legend
nva.eff.int.data$type <- factor(nva.eff.int.data$type, levels = c("0", "1"))
nva.eff.int.data$party_nva <- factor(nva.eff.int.data$party_nva, levels = c("0", "1"))

# Save the interaction plot to a high-resolution PNG file
png("nva_correct_int_g.png", units="in", width=6, height=4, res=1200)

# Create a ggplot showing predicted probabilities of correct forecast by type and party preference
# --- Plotting: figure creation starts here. Review aesthetics (aes) and any facet/scales for readability.
nva.eff.int.g <- ggplot(nva.eff.int.data, aes(x=type, y=fit, color=party_nva, group=party_nva)) +
  geom_line(size=.5, show.legend = TRUE, aes(linetype=party_nva, color=party_nva)) +      # Lines for each party group
  geom_point(size=.5, show.legend = FALSE) +                                             # Points for predicted values
  geom_errorbar(aes(ymin=fit-se, ymax=fit+se, color=party_nva), width = 0.05, show.legend = FALSE) +  # Confidence intervals
  labs(title = "N-VA") +                                                                # Plot title
  scale_x_discrete(labels = c("0" = "Citizens", "1" = "Politicians")) +                  # Label x-axis categories
  scale_y_continuous(name = "Predicted Probability", limits=c(0, 1),                      # y-axis scale and label
                     breaks = seq(0, 1, .2)) +
  facet_wrap(~ survey, ncol = 2) +
  theme_minimal() +                                                                      # Minimal theme for clean look
  theme(axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9),
        axis.title = element_text(size = 9), 
        axis.title.y = element_text(margin = margin(t = 0, r = 5, b = 0, l = 0)),
        axis.title.x = element_blank(),
        legend.title = element_text(size = 8), legend.text=element_text(size=8),
        legend.key.size = unit(.8, 'cm'),
        plot.title = element_blank(),
        plot.margin=unit(c(0,.01,.3,0), "cm"), legend.position="bottom") +
  scale_colour_manual(values=c("#636363", "#FCBD1B"), name="Party Preference",            # Custom colors for party groups
                      labels=c("0"="Other","1"="N-VA"), guide = guide_legend(reverse=TRUE)) +
  scale_fill_manual(values=c("#636363", "#FCBD1B"), name="Party Preference",              # Fill colors (if needed)
                    labels=c("0"="Other","1"="N-VA"), guide = guide_legend(reverse=TRUE)) +
  scale_linetype_manual(values=c("longdash","solid"), name="Party Preference",            # Line types for party groups
                        labels=c("0"="Other","1"="N-VA"), guide = guide_legend(reverse=TRUE))

# Print the plot object to the file
nva.eff.int.g

# Close the PNG device to save the file
dev.off()

# =============================================================================
# 21. Interaction: Correct seat forecasts (Vlaams Belang)
# =============================================================================

# Fit logistic regression model
model.seats_vb_correct <- glm(seats_vb_correct ~ party_vb*type*survey + 
                                sex + age, family="binomial", data = flanders.data)

# Display model summary
summary(model.seats_vb_correct)

# Predicted probabilities of correct forecast by type and party preference
# Not computed: Nearly VB respondents predicted more seats for VB

# Predicted probabilities of correct forecast by party preference only
vb.pred.party.prop <- ggeffect(model.seats_vb_correct, terms = c("party_vb"))
vb.pred.party.ref <- ggpredict(model.seats_vb_correct, terms = c("party_vb"))

# Predicted probabilities of correct forecast by respondent type only (citizen or politician)
vb.pred.type.prop <- ggeffect(model.seats_vb_correct, terms = c("type"))
vb.pred.type.ref <- ggpredict(model.seats_vb_correct, terms = c("type"))

# Compute Average Marginal Effects (AMEs) for the interaction between party preference and type
vb.eff.int <- effect(term="party_vb*type*survey", mod=model.seats_vb_correct)

# Convert the effect object to a data frame for easier manipulation and plotting
vb.eff.int.data <- as.data.frame(vb.eff.int)

# Add a column labeling the party for clarity in plots/tables
vb.eff.int.data$party <- "VB (+8)"

# Show the effect data frame (predicted probabilities and confidence intervals)
vb.eff.int.data

# Test significance of the interaction between party preference and type using estimated marginal means (EMMs)
vb.int.f <- emmeans(model.seats_vb_correct, ~ type | party_vb * survey)
vb.int.s <- pairs(pairs(emmeans::regrid(vb.int.f)), by = NULL)

# Plot the first difference effects with comparison indicators
plot(vb.int.f, comparisons = TRUE)

# Prepare factor levels for plotting to ensure correct order on x-axis and legend
vb.eff.int.data$type <- factor(vb.eff.int.data$type, levels = c("0", "1"))
vb.eff.int.data$party_vb <- factor(vb.eff.int.data$party_vb, levels = c("0", "1"))

# Save the interaction plot to a high-resolution PNG file
png("vb_correct_int_g.png", units="in", width=6, height=4, res=1200)

# Create a ggplot showing predicted probabilities of correct forecast by type and party preference
# --- Plotting: figure creation starts here. Review aesthetics (aes) and any facet/scales for readability.
vb.eff.int.g <- ggplot(vb.eff.int.data, aes(x=type, y=fit, color=party_vb, group=party_vb)) +
  geom_line(size=.5, show.legend = TRUE, aes(linetype=party_vb, color=party_vb)) +      # Lines for each party group
  geom_point(size=.5, show.legend = FALSE) +                                             # Points for predicted values
  geom_errorbar(aes(ymin=fit-se, ymax=fit+se, color=party_vb), width = 0.05, show.legend = FALSE) +  # Confidence intervals
  labs(title = "Vlaams Belang") +                                                                # Plot title
  scale_x_discrete(labels = c("0" = "Citizens", "1" = "Politicians")) +                  # Label x-axis categories
  scale_y_continuous(name = "Predicted Probability", limits=c(0, 1),                      # y-axis scale and label
                     breaks = seq(0, 1, .2)) +
  facet_wrap(~ survey, ncol = 2) +
  theme_minimal() +                                                                      # Minimal theme for clean look
  theme(axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9),
        axis.title = element_text(size = 9), 
        axis.title.y = element_text(margin = margin(t = 0, r = 5, b = 0, l = 0)),
        axis.title.x = element_blank(),
        legend.title = element_text(size = 8), legend.text=element_text(size=8),
        legend.key.size = unit(.8, 'cm'),
        plot.title = element_blank(),
        plot.margin=unit(c(0,.01,.3,0), "cm"), legend.position="bottom") +
  scale_colour_manual(values=c("#636363", "black"), name="Party Preference",            # Custom colors for party groups
                      labels=c("0"="Other","1"="Vlaams Belang"), guide = guide_legend(reverse=TRUE)) +
  scale_fill_manual(values=c("#636363", "black"), name="Party Preference",              # Fill colors (if needed)
                    labels=c("0"="Other","1"="Vlaams Belang"), guide = guide_legend(reverse=TRUE)) +
  scale_linetype_manual(values=c("longdash","solid"), name="Party Preference",            # Line types for party groups
                        labels=c("0"="Other","1"="Vlaams Belang"), guide = guide_legend(reverse=TRUE))

# Print the plot object to the file
vb.eff.int.g

# Close the PNG device to save the file
dev.off()

# =============================================================================
# 22. Interaction: Correct seat forecasts - Merge all AMES
# =============================================================================

# Standardize column names: rename party-specific binary indicator columns to 'inparty'
# This simplifies merging and plotting across different parties
pvda.eff.int.data   <- pvda.eff.int.data   %>% dplyr::rename(inparty = party_pvda)
groen.eff.int.data   <- groen.eff.int.data   %>% dplyr::rename(inparty = party_groen)
vooruit.eff.int.data <- vooruit.eff.int.data %>% dplyr::rename(inparty = party_vooruit)
cdv.eff.int.data     <- cdv.eff.int.data     %>% dplyr::rename(inparty = party_cdv)
ovld.eff.int.data    <- ovld.eff.int.data    %>% dplyr::rename(inparty = party_ovld)
nva.eff.int.data     <- nva.eff.int.data     %>% dplyr::rename(inparty = party_nva)
vb.eff.int.data   <- vb.eff.int.data   %>% dplyr::rename(inparty = party_vb)

# Combine (stack) the data frames into a single data frame for plotting
eff.int.data <- rbind(
  pvda.eff.int.data,
  groen.eff.int.data,
  vooruit.eff.int.data,
  cdv.eff.int.data,
  ovld.eff.int.data,
  nva.eff.int.data,
  vb.eff.int.data
)

# Reorder party levels for plotting (manual order based on result narrative or effect sizes)
eff.int.data$party <- factor(eff.int.data$party, levels = c(
  "PVDA (+5)",
  "Open Vld (–7)", 
  "Groen (–5)", 
  "N-VA (–4)", 
  "CD&V (–3)",
  "Vooruit (+5)",
  "VB (+8)"
))

# =============================================================================
# 23. Correct seat forecasts: Plot party preference x politician interaction
# =============================================================================

# Create visualization: Interaction effects with triple interaction
png("eff_int_g.png", units="in", width=5, height=7, res=1200)

# --- Plotting: figure creation starts here. Review aesthetics (aes) and any facet/scales for readability.
eff.int.g <- ggplot(eff.int.data, aes(x=type, y=fit, color=inparty, group=inparty)) +
  geom_point(size = .5, show.legend = FALSE) +
  geom_line(size = .5, aes(linetype = inparty)) +
  geom_ribbon(aes(ymin = fit - se, ymax = fit + se, fill = inparty),
              alpha = 0.3, linetype = 0, show.legend = FALSE) +
  
  scale_x_discrete(labels = c("0" = "Cit.", "1" = "Pol."),
                   expand = expansion(add = c(.1, .1))) +
  
  scale_y_continuous(name = "Predicted Probability", limits = c(0, 1),
                     breaks = seq(0, 1, .2),
                     expand = expansion(add = c(.1, .1))) +
  
  theme_bw() +
  theme(
    strip.text = element_text(size = 8),   # facet titles (both rows + columns)
    
    axis.text.x        = element_text(size = 9),
    axis.text.y        = element_text(size = 8),
    axis.title         = element_text(size = 10),
    axis.title.y       = element_text(size = 10, margin = margin(t = 0, r = 8, b = 0, l = 0)),
    axis.title.x       = element_blank(),
    
    legend.title       = element_text(size = 8),
    legend.text        = element_text(size = 8),
    legend.key.size    = unit(.6, "cm"),
    
    legend.position    = "bottom",
    legend.box         = "horizontal",
    legend.direction   = "horizontal",
    legend.spacing.x   = unit(0.6, "cm"),
    
    plot.title         = element_text(size = 12, hjust = 0.5, colour = "black"),
    panel.grid.major   = element_blank(),
    panel.grid.minor   = element_blank(),
    plot.margin        = unit(c(0, .01, .01, 0), "cm")
  ) +
  
  scale_colour_manual(
    values = c("#F8766D", "#00BFC4"),
    name = "Party Supporter?",
    labels = c("No", "Yes"),
    guide = guide_legend(title.position = "left", nrow = 1, byrow = TRUE)
  ) +
  scale_fill_manual(
    values = c("#F8766D", "#00BFC4"),
    name = "Party Supporter?",
    labels = c("No", "Yes"),
    guide = "none"
  ) +
  scale_linetype_manual(
    values = c("longdash", "solid"),
    name = "Party Supporter?",
    labels = c("No", "Yes"),
    guide = guide_legend(title.position = "top", nrow = 1, byrow = TRUE)
  ) +
  
  facet_grid(party ~ survey)

eff.int.g
dev.off()

# =============================================================================
# 24. Linear prediction: Probability of being in government (PVDA)
# =============================================================================

# Create a new data frame with all combinations of key predictors for prediction
# - party_pvda: binary indicator for PVDA supporter ("0" = no, "1" = yes)
# - type: respondent type ("0" = citizen, "1" = politician)
# Other predictors (ideology, age, sex, survey) are set to their sample means or modal category to hold constant
model.pvda.df <- expand.grid(
  party_pvda = c("0", "1"),                          # Party support indicator levels
  type = c("0", "1"),                               # Respondent type levels
  ideology = mean(flanders.data$ideology, na.rm=TRUE), # Continuous covariate set to mean
  age = mean(flanders.data$age, na.rm=TRUE),           # Continuous covariate set to mean
  sex = modal(flanders.data$sex, na.rm=TRUE),           # Categorical covariate set to mode
  survey = unique(flanders.data$survey)      # Both survey years
)

# Use the previously fitted model 'model.govt_pvda' to predict probability of being in government
# type = "response" returns predicted probabilities on the original scale (0 to 1)
# se.fit = TRUE returns standard errors of the predictions
# interval = "confidence" and level = 0.95 for 95% confidence intervals (may or may not be used below)
model.govt.pvda.pr <- predict(
  model.govt_pvda, 
  newdata = model.pvda.df, 
  type = "response", 
  se.fit = TRUE, 
  interval = "confidence",
  level = 0.95
)

# Extract the predicted fit values and standard errors into a data frame for easier manipulation
model.govt.pvda.pr.fit <- as.data.frame(model.govt.pvda.pr$fit)

# Add a row identifier to both data frames to enable merging
model.pvda.df <- model.pvda.df %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  dplyr::mutate(row_id = row_number())

model.govt.pvda.pr.fit <- model.govt.pvda.pr.fit %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  dplyr::mutate(row_id = row_number())

# Merge the predictions back with the original prediction grid data by row_id
model.govt.pvda.df <- merge(model.pvda.df, model.govt.pvda.pr.fit, by = "row_id")

# Add a new column labeling the model/source for clarity in combined datasets or plotting
model.govt.pvda.df$model <- "PVDA"

# Rename the party_pvda column to a more generic 'party' for consistency with other datasets
model.govt.pvda.df <- model.govt.pvda.df %>%
  dplyr::rename(party = party_pvda)

# Open a PNG device to save the plot as a high-resolution image
png("model_govt_pvda_g.png", units="in", width=12, height=8, res=1200)

# Create the ggplot object using the model predictions dataframe
# --- Plotting: figure creation starts here. Review aesthetics (aes) and any facet/scales for readability.
model.govt.pvda.g <- ggplot(model.govt.pvda.df, aes(x = factor(type), y = fit, 
                                                    color = factor(party), group = party, shape = party)) +
  # Add lines for each party (grouped by party), with line type mapped to party
  geom_line(aes(linetype = party)) +
  # Add points for each prediction; shape mapped to party
  geom_point(aes(shape = factor(party)), 
             # position = position_dodge(0.1),  # Optionally adjust horizontal position
             size = 3.5) +  # Size of points
  # Add error bars for confidence intervals (lower = lwr, upper = upr)
  geom_errorbar(aes(ymin = lwr, ymax = upr), 
                # position = position_dodge(0.1),  # Optionally dodge error bars
                width = 0.05) +  # Narrow width for error bars
  # Customize the x-axis: categorical with labels for 0 and 1
  scale_x_discrete(name = "", breaks = c(0, 1), labels = c("Citizens", "Politicians")) +
  # Label for the y-axis
  scale_y_continuous(name = "Linear Prediction") +
  facet_wrap(~ survey, ncol = 2) +
  # Use a minimal theme for a clean plot look
  theme_minimal() +
  # Customize text and layout of axes and legend
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = 1),
    panel.spacing.x = unit(1.5, "lines"),  # Space between facets
    axis.title.x = element_text(size = 20),  # X-axis title font size
    axis.title.y = element_text(size = 20, margin = margin(r = 12)),  # Y-axis title with right margin
    axis.text.x = element_text(size = 18),   # X-axis tick labels
    axis.text.y = element_text(size = 18),   # Y-axis tick labels
    strip.text = element_text(size = 18),
    legend.position = "bottom",         
    legend.text = element_text(size = 17),   # Legend text size
    legend.title = element_text(size = 18),   # Legend title size
    legend.key.height = unit(.8, "lines"),   # Legend key height
    legend.key.width = unit(1, 'cm')         # Legend key width
  ) +
  # Manually define colors for each party
  scale_color_manual("Party Preference", values = c("#636363", "#990000"), labels = party.labels.pvda) +
  # Manually define line types for each party
  scale_linetype_manual("Party Preference", values = c("dashed", "solid"), labels = party.labels.pvda) +
  # Manually define point shapes for each party
  scale_shape_manual("Party Preference", values = c(16, 17), labels = party.labels.pvda)

# Print the plot to the PNG device
model.govt.pvda.g

# Close the graphics device and save the file
dev.off()

# =============================================================================
# 25. Linear prediction: Probability OF being government (Groen)
# =============================================================================

# Create a new data frame with all combinations of key predictors for prediction
# - party_groen: binary indicator for Groen supporter ("0" = no, "1" = yes)
# - type: respondent type ("0" = citizen, "1" = politician)
# Other predictors (ideology, age, sex, survey) are set to their sample means or modal category to hold constant
model.groen.df <- expand.grid(
  party_groen = c("0", "1"),
  type = c("0", "1"),
  ideology = mean(flanders.data$ideology, na.rm=TRUE),
  age = mean(flanders.data$age, na.rm=TRUE),
  sex = modal(flanders.data$sex, na.rm=TRUE),
  survey = unique(flanders.data$survey)
)

# Use the previously fitted model 'model.govt_groen' to predict probability of being in government
# type = "response" returns predicted probabilities on the original scale (0 to 1)
# se.fit = TRUE returns standard errors of the predictions
# interval = "confidence" and level = 0.95 for 95% confidence intervals
model.govt.groen.pr <- predict(
  model.govt_groen, 
  newdata = model.groen.df, 
  type = "response", 
  se.fit = TRUE, 
  interval = "confidence",
  level = 0.95
)

# Extract the predicted fit values and standard errors into a data frame for easier manipulation
model.govt.groen.pr.fit <- as.data.frame(model.govt.groen.pr$fit)

# Add a row identifier to both data frames to enable merging
model.groen.df <- model.groen.df %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  dplyr::mutate(row_id = row_number())

model.govt.groen.pr.fit <- model.govt.groen.pr.fit %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  dplyr::mutate(row_id = row_number())

# Merge the predictions back with the original prediction grid data by row_id
model.govt.groen.df <- merge(model.groen.df, model.govt.groen.pr.fit, by = "row_id")

# Add a new column labeling the model/source for clarity in combined datasets or plotting
model.govt.groen.df$model <- "Groen"

# Rename the party_groen column to a more generic 'party' for consistency with other datasets
model.govt.groen.df <- model.govt.groen.df %>%
  dplyr::rename(party = party_groen)

# Open a PNG device to save the plot as a high-resolution image
png("model_govt_groen_g.png", units="in", width=12, height=8, res=1200)

# Create the ggplot object using the model predictions dataframe
# --- Plotting: figure creation starts here. Review aesthetics (aes) and any facet/scales for readability.
model.govt.groen.g <- ggplot(model.govt.groen.df, aes(x = factor(type), y = fit, 
                                                      color = factor(party), group = party, shape = party)) +
  # Add lines for each party (grouped by party), with line type mapped to party
  geom_line(aes(linetype = party)) +
  # Add points for each prediction; shape mapped to party
  geom_point(aes(shape = factor(party)), 
             size = 3.5) +
  # Add error bars for confidence intervals (lower = lwr, upper = upr)
  geom_errorbar(aes(ymin = lwr, ymax = upr), 
                width = 0.05) +
  # Customize the x-axis: categorical with labels for 0 and 1
  scale_x_discrete(name = "", breaks = c(0, 1), labels = c("Citizens", "Politicians")) +
  # Label for the y-axis
  scale_y_continuous(name = "Linear Prediction") +
  # Use a minimal theme for a clean plot look
  facet_wrap(~ survey, ncol = 2) +
  theme_minimal() +
  # Customize text and layout of axes and legend
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = 1),
    panel.spacing.x = unit(1.5, "lines"),  # Space between facets
    axis.title.x = element_text(size = 20),
    axis.title.y = element_text(size = 20, margin = margin(r = 12)),
    axis.text.x = element_text(size = 18),
    axis.text.y = element_text(size = 18),
    strip.text = element_text(size = 18),
    legend.position = "bottom",         
    legend.text = element_text(size = 17),   # Legend text size
    legend.title = element_text(size = 18),   # Legend title size
    legend.key.height = unit(.8, "lines"),
    legend.key.width = unit(1, 'cm')
  ) +
  # Manually define colors for each party
  scale_color_manual("Party Preference", values = c("#636363", "#04847A"), labels = party.labels.groen) +
  # Manually define line types for each party
  scale_linetype_manual("Party Preference", values = c("dashed", "solid"), labels = party.labels.groen) +
  # Manually define point shapes for each party
  scale_shape_manual("Party Preference", values = c(16, 17), labels = party.labels.groen)

# Print the plot to the PNG device
model.govt.groen.g

# Close the graphics device and save the file
dev.off()

# =============================================================================
# 26. Linear prediction: Probability of being in government (Vooruit)
# =============================================================================

# Create a new data frame with all combinations of key predictors for prediction
# - party_vooruit: binary indicator for Vooruit supporter ("0" = no, "1" = yes)
# - type: respondent type ("0" = citizen, "1" = politician)
# Other predictors (ideology, age, sex, survey) are set to their sample means or modal category to hold constant
model.vooruit.df <- expand.grid(
  party_vooruit = c("0", "1"),
  type = c("0", "1"),
  ideology = mean(flanders.data$ideology, na.rm=TRUE),
  age = mean(flanders.data$age, na.rm=TRUE),
  sex = modal(flanders.data$sex, na.rm=TRUE),
  survey = unique(flanders.data$survey)
)

# Use the previously fitted model 'model.govt_vooruit' to predict probability of being in government
# type = "response" returns predicted probabilities on the original scale (0 to 1)
# se.fit = TRUE returns standard errors of the predictions
# interval = "confidence" and level = 0.95 for 95% confidence intervals
model.govt.vooruit.pr <- predict(
  model.govt_vooruit, 
  newdata = model.vooruit.df, 
  type = "response", 
  se.fit = TRUE, 
  interval = "confidence",
  level = 0.95
)

# Extract the predicted fit values and standard errors into a data frame for easier manipulation
model.govt.vooruit.pr.fit <- as.data.frame(model.govt.vooruit.pr$fit)

# Add a row identifier to both data frames to enable merging
model.vooruit.df <- model.vooruit.df %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  dplyr::mutate(row_id = row_number())

model.govt.vooruit.pr.fit <- model.govt.vooruit.pr.fit %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  dplyr::mutate(row_id = row_number())

# Merge the predictions back with the original prediction grid data by row_id
model.govt.vooruit.df <- merge(model.vooruit.df, model.govt.vooruit.pr.fit, by = "row_id")

# Add a new column labeling the model/source for clarity in combined datasets or plotting
model.govt.vooruit.df$model <- "Vooruit"

# Rename the party_vooruit column to a more generic 'party' for consistency with other datasets
model.govt.vooruit.df <- model.govt.vooruit.df %>%
  dplyr::rename(party = party_vooruit)

# Open a PNG device to save the plot as a high-resolution image
png("model_govt_vooruit_g.png", units="in", width=12, height=8, res=1200)

# Create the ggplot object using the model predictions dataframe
# --- Plotting: figure creation starts here. Review aesthetics (aes) and any facet/scales for readability.
model.govt.vooruit.g <- ggplot(model.govt.vooruit.df, aes(x = factor(type), y = fit, 
                                                          color = factor(party), group = party, shape = party)) +
  # Add lines for each party (grouped by party), with line type mapped to party
  geom_line(aes(linetype = party)) +
  # Add points for each prediction; shape mapped to party
  geom_point(aes(shape = factor(party)), 
             size = 3.5) +
  # Add error bars for confidence intervals (lower = lwr, upper = upr)
  geom_errorbar(aes(ymin = lwr, ymax = upr), 
                width = 0.05) +
  # Customize the x-axis: categorical with labels for 0 and 1
  scale_x_discrete(name = "", breaks = c(0, 1), labels = c("Citizens", "Politicians")) +
  # Label for the y-axis
  scale_y_continuous(name = "Linear Prediction") +
  # Use a minimal theme for a clean plot look
  facet_wrap(~ survey, ncol = 2) +
  theme_minimal() +
  # Customize text and layout of axes and legend
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = 1),
    panel.spacing.x = unit(1.5, "lines"),  # Space between facets
    axis.title.x = element_text(size = 20),
    axis.title.y = element_text(size = 20, margin = margin(r = 12)),
    axis.text.x = element_text(size = 18),
    axis.text.y = element_text(size = 18),
    strip.text = element_text(size = 18),
    legend.position = "bottom",         
    legend.text = element_text(size = 17),   # Legend text size
    legend.title = element_text(size = 18),   # Legend title size
    legend.key.height = unit(.8, "lines"),
    legend.key.width = unit(1, 'cm')
  ) +
  # Manually define colors for each party
  scale_color_manual("Party Preference", values = c("#636363", "#FF0000"), labels = party.labels.vooruit) +
  # Manually define line types for each party
  scale_linetype_manual("Party Preference", values = c("dashed", "solid"), labels = party.labels.vooruit) +
  # Manually define point shapes for each party
  scale_shape_manual("Party Preference", values = c(16, 17), labels = party.labels.vooruit)

# Print the plot to the PNG device
model.govt.vooruit.g

# Close the graphics device and save the file
dev.off()

# =============================================================================
# 27. Linear prediction: Probability of being in government (CD&V)
# =============================================================================

# Create a new data frame with all combinations of key predictors for prediction
# - party_cdv: binary indicator for CD&V supporter ("0" = no, "1" = yes)
# - type: respondent type ("0" = citizen, "1" = politician)
# Other predictors (ideology, age, sex, survey) are set to their sample means or modal category to hold constant
model.cdv.df <- expand.grid(
  party_cdv = c("0", "1"),
  type = c("0", "1"),
  ideology = mean(flanders.data$ideology, na.rm=TRUE),
  age = mean(flanders.data$age, na.rm=TRUE),
  sex = modal(flanders.data$sex, na.rm=TRUE),
  survey = unique(flanders.data$survey)
)

# Use the previously fitted model 'model.govt_cdv' to predict probability of being in government
# type = "response" returns predicted probabilities on the original scale (0 to 1)
# se.fit = TRUE returns standard errors of the predictions
# interval = "confidence" and level = 0.95 for 95% confidence intervals
model.govt.cdv.pr <- predict(
  model.govt_cdv, 
  newdata = model.cdv.df, 
  type = "response", 
  se.fit = TRUE, 
  interval = "confidence",
  level = 0.95
)

# Extract the predicted fit values and standard errors into a data frame for easier manipulation
model.govt.cdv.pr.fit <- as.data.frame(model.govt.cdv.pr$fit)

# Add a row identifier to both data frames to enable merging
model.cdv.df <- model.cdv.df %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  dplyr::mutate(row_id = row_number())

model.govt.cdv.pr.fit <- model.govt.cdv.pr.fit %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  dplyr::mutate(row_id = row_number())

# Merge the predictions back with the original prediction grid data by row_id
model.govt.cdv.df <- merge(model.cdv.df, model.govt.cdv.pr.fit, by = "row_id")

# Add a new column labeling the model/source for clarity in combined datasets or plotting
model.govt.cdv.df$model <- "CD&V"

# Rename the party_cdv column to a more generic 'party' for consistency with other datasets
model.govt.cdv.df <- model.govt.cdv.df %>%
  dplyr::rename(party = party_cdv)

# Open a PNG device to save the plot as a high-resolution image
png("model_govt_cdv_g.png", units="in", width=12, height=8, res=1200)

# Create the ggplot object using the model predictions dataframe
# --- Plotting: figure creation starts here. Review aesthetics (aes) and any facet/scales for readability.
model.govt.cdv.g <- ggplot(model.govt.cdv.df, aes(x = factor(type), y = fit, 
                                                  color = factor(party), group = party, shape = party)) +
  # Add lines for each party (grouped by party), with line type mapped to party
  geom_line(aes(linetype = party)) +
  # Add points for each prediction; shape mapped to party
  geom_point(aes(shape = factor(party)), 
             size = 3.5) +
  # Add error bars for confidence intervals (lower = lwr, upper = upr)
  geom_errorbar(aes(ymin = lwr, ymax = upr), 
                width = 0.05) +
  # Customize the x-axis: categorical with labels for 0 and 1
  scale_x_discrete(name = "", breaks = c(0, 1), labels = c("Citizens", "Politicians")) +
  # Label for the y-axis
  scale_y_continuous(name = "Linear Prediction") +
  facet_wrap(~ survey, ncol = 2) +
  # Use a minimal theme for a clean plot look
  theme_minimal() +
  # Customize text and layout of axes and legend
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = 1),
    panel.spacing.x = unit(1.5, "lines"),  # Space between facets
    axis.title.x = element_text(size = 20),
    axis.title.y = element_text(size = 20, margin = margin(r = 12)),
    axis.text.x = element_text(size = 18),
    axis.text.y = element_text(size = 18),
    strip.text = element_text(size = 18),
    legend.position = "bottom",         
    legend.text = element_text(size = 17),   # Legend text size
    legend.title = element_text(size = 18),   # Legend title size
    legend.key.height = unit(.8, "lines"),
    legend.key.width = unit(1, 'cm')
  ) +
  # Manually define colors for each party
  scale_color_manual("Party Preference", values = c("#636363", "#FF6000"), labels = party.labels.cdv) +
  # Manually define line types for each party
  scale_linetype_manual("Party Preference", values = c("dashed", "solid"), labels = party.labels.cdv) +
  # Manually define point shapes for each party
  scale_shape_manual("Party Preference", values = c(16, 17), labels = party.labels.cdv)

# Print the plot to the PNG device
model.govt.cdv.g

# Close the graphics device and save the file
dev.off()

# =============================================================================
# 28. Linear prediction: Probability of being in government (Open Vld)
# =============================================================================

# Create a new data frame with all combinations of key predictors for prediction
# - party_ovld: binary indicator for Open VLD supporter ("0" = no, "1" = yes)
# - type: respondent type ("0" = citizen, "1" = politician)
# Other predictors (ideology, age, sex, survey) are set to their sample means or modal category to hold constant
model.ovld.df <- expand.grid(
  party_ovld = c("0", "1"),
  type = c("0", "1"),
  ideology = mean(flanders.data$ideology, na.rm=TRUE),
  age = mean(flanders.data$age, na.rm=TRUE),
  sex = modal(flanders.data$sex, na.rm=TRUE),
  survey = unique(flanders.data$survey)
)

# Use the previously fitted model 'model.govt_ovld' to predict probability of being in government
# type = "response" returns predicted probabilities on the original scale (0 to 1)
# se.fit = TRUE returns standard errors of the predictions
# interval = "confidence" and level = 0.95 for 95% confidence intervals
model.govt.ovld.pr <- predict(
  model.govt_ovld, 
  newdata = model.ovld.df, 
  type = "response", 
  se.fit = TRUE, 
  interval = "confidence",
  level = 0.95
)

# Extract the predicted fit values and standard errors into a data frame for easier manipulation
model.govt.ovld.pr.fit <- as.data.frame(model.govt.ovld.pr$fit)

# Add a row identifier to both data frames to enable merging
model.ovld.df <- model.ovld.df %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  dplyr::mutate(row_id = row_number())

model.govt.ovld.pr.fit <- model.govt.ovld.pr.fit %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  dplyr::mutate(row_id = row_number())

# Merge the predictions back with the original prediction grid data by row_id
model.govt.ovld.df <- merge(model.ovld.df, model.govt.ovld.pr.fit, by = "row_id")

# Add a new column labeling the model/source for clarity in combined datasets or plotting
model.govt.ovld.df$model <- "Open Vld"

# Rename the party_ovld column to a more generic 'party' for consistency with other datasets
model.govt.ovld.df <- model.govt.ovld.df %>%
  dplyr::rename(party = party_ovld)

# Open a PNG device to save the plot as a high-resolution image
png("model_govt_ovld_g.png", units="in", width=12, height=8, res=1200)

# Create the ggplot object using the model predictions dataframe
# --- Plotting: figure creation starts here. Review aesthetics (aes) and any facet/scales for readability.
model.govt.ovld.g <- ggplot(model.govt.ovld.df, aes(x = factor(type), y = fit, 
                                                    color = factor(party), group = party, shape = party)) +
  # Add lines for each party (grouped by party), with line type mapped to party
  geom_line(aes(linetype = party)) +
  # Add points for each prediction; shape mapped to party
  geom_point(aes(shape = factor(party)), 
             size = 3.5) +
  # Add error bars for confidence intervals (lower = lwr, upper = upr)
  geom_errorbar(aes(ymin = lwr, ymax = upr), 
                width = 0.05) +
  # Customize the x-axis: categorical with labels for 0 and 1
  scale_x_discrete(name = "", breaks = c(0, 1), labels = c("Citizens", "Politicians")) +
  # Label for the y-axis
  scale_y_continuous(name = "Linear Prediction") +
  facet_wrap(~ survey, ncol = 2) +
  # Use a minimal theme for a clean plot look
  theme_minimal() +
  # Customize text and layout of axes and legend
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = 1),
    panel.spacing.x = unit(1.5, "lines"),  # Space between facets
    axis.title.x = element_text(size = 20),
    axis.title.y = element_text(size = 20, margin = margin(r = 12)),
    axis.text.x = element_text(size = 18),
    axis.text.y = element_text(size = 18),
    strip.text = element_text(size = 18),
    legend.position = "bottom",         
    legend.text = element_text(size = 17),   # Legend text size
    legend.title = element_text(size = 18),   # Legend title size
    legend.key.height = unit(.8, "lines"),
    legend.key.width = unit(1, 'cm')
  ) +
  # Manually define colors for each party
  scale_color_manual("Party Preference", values = c("#636363", "#005DAA"), labels = party.labels.ovld) +
  # Manually define line types for each party
  scale_linetype_manual("Party Preference", values = c("dashed", "solid"), labels = party.labels.ovld) +
  # Manually define point shapes for each party
  scale_shape_manual("Party Preference", values = c(16, 17), labels = party.labels.ovld)

# Print the plot to the PNG device
model.govt.ovld.g

# Close the graphics device and save the file
dev.off()

# =============================================================================
# 29. Linear prediction: Probability of being in government (N-VA)
# =============================================================================

# Create a new data frame with all combinations of key predictors for prediction
# - party_nva: binary indicator for N-VA supporter ("0" = no, "1" = yes)
# - type: respondent type ("0" = citizen, "1" = politician)
# Other predictors (ideology, age, sex, survey) are set to their sample means or modal category to hold constant
model.nva.df <- expand.grid(
  party_nva = c("0", "1"),
  type = c("0", "1"),
  ideology = mean(flanders.data$ideology, na.rm=TRUE),
  age = mean(flanders.data$age, na.rm=TRUE),
  sex = modal(flanders.data$sex, na.rm=TRUE),
  survey = unique(flanders.data$survey)
)

# Use the previously fitted model 'model.govt_nva' to predict probability of being in government
# type = "response" returns predicted probabilities on the original scale (0 to 1)
# se.fit = TRUE returns standard errors of the predictions
# interval = "confidence" and level = 0.95 for 95% confidence intervals
model.govt.nva.pr <- predict(
  model.govt_nva, 
  newdata = model.nva.df, 
  type = "response", 
  se.fit = TRUE, 
  interval = "confidence",
  level = 0.95
)

# Extract the predicted fit values and standard errors into a data frame for easier manipulation
model.govt.nva.pr.fit <- as.data.frame(model.govt.nva.pr$fit)

# Add a row identifier to both data frames to enable merging
model.nva.df <- model.nva.df %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  dplyr::mutate(row_id = row_number())

model.govt.nva.pr.fit <- model.govt.nva.pr.fit %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  dplyr::mutate(row_id = row_number())

# Merge the predictions back with the original prediction grid data by row_id
model.govt.nva.df <- merge(model.nva.df, model.govt.nva.pr.fit, by = "row_id")

# Add a new column labeling the model/source for clarity in combined datasets or plotting
model.govt.nva.df$model <- "N-VA"

# Rename the party_nva column to a more generic 'party' for consistency with other datasets
model.govt.nva.df <- model.govt.nva.df %>%
  dplyr::rename(party = party_nva)

# Open a PNG device to save the plot as a high-resolution image
png("model_govt_nva_g.png", units="in", width=12, height=8, res=1200)

# Create the ggplot object using the model predictions dataframe
# --- Plotting: figure creation starts here. Review aesthetics (aes) and any facet/scales for readability.
model.govt.nva.g <- ggplot(model.govt.nva.df, aes(x = factor(type), y = fit, 
                                                  color = factor(party), group = party, shape = party)) +
  # Add lines for each party (grouped by party), with line type mapped to party
  geom_line(aes(linetype = party)) +
  # Add points for each prediction; shape mapped to party
  geom_point(aes(shape = factor(party)), 
             size = 3.5) +
  # Add error bars for confidence intervals (lower = lwr, upper = upr)
  geom_errorbar(aes(ymin = lwr, ymax = upr), 
                width = 0.05) +
  # Customize the x-axis: categorical with labels for 0 and 1
  scale_x_discrete(name = "", breaks = c(0, 1), labels = c("Citizens", "Politicians")) +
  # Label for the y-axis
  scale_y_continuous(name = "Linear Prediction") +
  facet_wrap(~ survey, ncol = 2) +
  # Use a minimal theme for a clean plot look
  theme_minimal() +
  # Customize text and layout of axes and legend
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = 1),
    panel.spacing.x = unit(1.5, "lines"),  # Space between facets
    axis.title.x = element_text(size = 20),
    axis.title.y = element_text(size = 20, margin = margin(r = 12)),
    axis.text.x = element_text(size = 18),
    axis.text.y = element_text(size = 18),
    strip.text = element_text(size = 18),
    legend.position = "bottom",         
    legend.text = element_text(size = 17),   # Legend text size
    legend.title = element_text(size = 18),   # Legend title size
    legend.key.height = unit(.8, "lines"),
    legend.key.width = unit(1, 'cm')
  ) +
  # Manually define colors for each party
  scale_color_manual("Party Preference", values = c("#636363", "#FCBD1B"), labels = party.labels.ovld) +
  # Manually define line types for each party
  scale_linetype_manual("Party Preference", values = c("dashed", "solid"), labels = party.labels.ovld) +
  # Manually define point shapes for each party
  scale_shape_manual("Party Preference", values = c(16, 17), labels = party.labels.ovld)

# Print the plot to the PNG device
model.govt.nva.g

# Close the graphics device and save the file
dev.off()

# =============================================================================
# 30. Linear prediction: Probability of being in government (Vlaams Belang)
# =============================================================================

# Create a new data frame with all combinations of key predictors for prediction
# - party_vb: binary indicator for Vlaams Belang supporter ("0" = no, "1" = yes)
# - type: respondent type ("0" = citizen, "1" = politician)
# Other predictors (ideology, age, sex, survey) are set to their sample means or modal category to hold constant
model.vb.df <- expand.grid(
  party_vb = c("0", "1"),
  type = c("0", "1"),
  ideology = mean(flanders.data$ideology, na.rm=TRUE),
  age = mean(flanders.data$age, na.rm=TRUE),
  sex = modal(flanders.data$sex, na.rm=TRUE),
  survey = unique(flanders.data$survey)
)

# Use the previously fitted model 'model.govt_vb' to predict probability of being in government
# type = "response" returns predicted probabilities on the original scale (0 to 1)
# se.fit = TRUE returns standard errors of the predictions
# interval = "confidence" and level = 0.95 for 95% confidence intervals
model.govt.vb.pr <- predict(
  model.govt_vb, 
  newdata = model.vb.df, 
  type = "response", 
  se.fit = TRUE, 
  interval = "confidence",
  level = 0.95
)

# Extract the predicted fit values and standard errors into a data frame for easier manipulation
model.govt.vb.pr.fit <- as.data.frame(model.govt.vb.pr$fit)

# Add a row identifier to both data frames to enable merging
model.vb.df <- model.vb.df %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  dplyr::mutate(row_id = row_number())

model.govt.vb.pr.fit <- model.govt.vb.pr.fit %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  dplyr::mutate(row_id = row_number())

# Merge the predictions back with the original prediction grid data by row_id
model.govt.vb.df <- merge(model.vb.df, model.govt.vb.pr.fit, by = "row_id")

# Add a new column labeling the model/source for clarity in combined datasets or plotting
model.govt.vb.df$model <- "Vlaams Belang"

# Rename the party_vb column to a more generic 'party' for consistency with other datasets
model.govt.vb.df <- model.govt.vb.df %>%
  dplyr::rename(party = party_vb)

# Open a PNG device to save the plot as a high-resolution image
png("model_govt_vb_g.png", units="in", width=12, height=8, res=1200)

# Create the ggplot object using the model predictions dataframe
# --- Plotting: figure creation starts here. Review aesthetics (aes) and any facet/scales for readability.
model.govt.vb.g <- ggplot(model.govt.vb.df, aes(x = factor(type), y = fit, 
                                                color = factor(party), group = party, shape = party)) +
  # Add lines for each party (grouped by party), with line type mapped to party
  geom_line(aes(linetype = party)) +
  # Add points for each prediction; shape mapped to party
  geom_point(aes(shape = factor(party)), 
             size = 3.5) +
  # Add error bars for confidence intervals (lower = lwr, upper = upr)
  geom_errorbar(aes(ymin = lwr, ymax = upr), 
                width = 0.05) +
  # Customize the x-axis: categorical with labels for 0 and 1
  scale_x_discrete(name = "", breaks = c(0, 1), labels = c("Citizens", "Politicians")) +
  # Label for the y-axis
  scale_y_continuous(name = "Linear Prediction") +
  facet_wrap(~ survey, ncol = 2) +
  # Use a minimal theme for a clean plot look
  theme_minimal() +
  # Customize text and layout of axes and legend
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = 1),
    panel.spacing.x = unit(1.5, "lines"),  # Space between facets
    axis.title.x = element_text(size = 20),
    axis.title.y = element_text(size = 20, margin = margin(r = 12)),
    axis.text.x = element_text(size = 18),
    axis.text.y = element_text(size = 18),
    strip.text = element_text(size = 18),
    legend.position = "bottom",         
    legend.text = element_text(size = 17),   # Legend text size
    legend.title = element_text(size = 18),   # Legend title size
    legend.key.height = unit(.8, "lines"),
    legend.key.width = unit(1, 'cm')
  ) +
  # Manually define colors for each party
  scale_color_manual("Party Preference", values = c("#636363", "black"), labels = party.labels.vb) +
  # Manually define line types for each party
  scale_linetype_manual("Party Preference", values = c("dashed", "solid"), labels = party.labels.vb) +
  # Manually define point shapes for each party
  scale_shape_manual("Party Preference", values = c(16, 17), labels = party.labels.vb)

# Print the plot to the PNG device
model.govt.vb.g

# Close the graphics device and save the file
dev.off()

# =============================================================================
# 31. Linear prediction: Merge
# =============================================================================

# Combine all individual party prediction data frames into one master dataframe
# This dataframe will contain predicted probabilities of being in government for all parties
model.govt.df <- rbind(
  model.govt.pvda.df,    # Predictions for PVDA
  model.govt.groen.df,   # Predictions for Groen
  model.govt.vooruit.df, # Predictions for Vooruit
  model.govt.ovld.df,    # Predictions for Open Vld
  model.govt.cdv.df,     # Predictions for CD&V
  model.govt.nva.df,     # Predictions for N-VA
  model.govt.vb.df       # Predictions for Vlaams Belang
)

# Open a PNG graphics device to save the plot to a high-resolution image file
png("model_govt_g.png", units="in", width=15, height=10, res=1200)

# Create a ggplot object to visualize the linear predictions by party and respondent type
# --- Plotting: figure creation starts here. Review aesthetics (aes) and any facet/scales for readability.
model.govt.g <- ggplot(model.govt.df, aes(
  x = factor(type),       # X-axis: respondent type (0 = Citizens, 1 = Politicians)
  y = fit,                # Y-axis: predicted linear prediction (fit)
  color = factor(party),  # Color points/lines by party (binary: supporter or not)
  group = party,          # Group lines by party
  shape = party,          # Point shape indicates party membership
  linetype = party        # Line type indicates party membership
)) +
  # Add lines connecting predictions for each party
  geom_line() +
  # Add points at prediction values with party-specific shapes
  geom_point(aes(shape = factor(party)), 
             #position = position_dodge(0.1),  # Optional: horizontally dodge points to avoid overlap
             size = 3.5) +  # Size of the points
  # Add vertical error bars representing confidence intervals (lower=lwr, upper=upr)
  geom_errorbar(aes(ymin = lwr, ymax = upr), 
                #position = position_dodge(0.1),  # Optional dodge to align with points
                width = 0.05) +  # Width of error bars
  # Define the x-axis as a discrete variable with labels
  scale_x_discrete(
    name = "", 
    breaks = c(0, 1), 
    labels = c("Citizens", "Politicians")
  ) +
  # Define y-axis scale limits and breaks, label axis
  scale_y_continuous(
    limits = c(0, 10),       # Y-axis range from 0 to 10
    breaks = seq(0, 10, 2), # Break ticks every 2 units
    name = "Linear Prediction"
  ) +
  # Facet plot by 'model' column (i.e., party), arranged in 2 rows,
  # with free y-axis scales per facet to accommodate differences in ranges
  facet_grid(survey ~ model) +
  # Use a clean minimal theme for better aesthetics
  theme_minimal() +
  # Customize text sizes and legend placement
  theme(
    axis.title.x = element_text(size = 20),            # X-axis title font size
    axis.title.y = element_text(size = 20, margin = margin(r = 12)), # Y-axis title with right margin
    axis.text.x = element_text(size = 18),             # X-axis tick label size
    axis.text.y = element_text(size = 18),             # Y-axis tick label size
    legend.position = c(0.85, 0.20),                    # Position legend near bottom-right corner
    legend.text = element_text(size = 16),             # Legend text size
    legend.title = element_blank(),                     # Remove legend title for clean look
    legend.key.height = unit(.8, "lines"),              # Legend key height adjustment
    legend.key.width = unit(1, 'cm'),                   # Legend key width adjustment
    strip.text = element_text(size = 18)                # Facet strip label font size
  ) +
  # Manually specify colors for party preference (adjust colors and labels as needed)
  scale_color_manual("Party Preference", values = c("#636363", "black"), labels = party.labels) +
  # Manually specify line types for party preference
  scale_linetype_manual("Party Preference", values = c("dashed", "solid"), labels = party.labels) +
  # Manually specify point shapes for party preference
  scale_shape_manual("Party Preference", values = c(16, 17), labels = party.labels)

# Print the plot to the PNG device (i.e., save the image)
model.govt.g

# Close the PNG device to finalize and save the image file
dev.off()

# =============================================================================
# 32. Predicted probabilities: Getting more or less seats (PVDA)
# =============================================================================

# Predict class probabilities using the trained model for PVDA
model.seats.pvda.pr <- predict(model.seats_pvda, newdata = model.pvda.df, 
                               type = "prob")  # Return probabilities for each class (fewer/same/more)

# Combine the original data with predicted probabilities
model.seats.pvda.int <- cbind(model.pvda.df, model.seats.pvda.pr)

# Reshape data from wide to long format for ggplot (e.g., fit.Fewer, fit.Same → rows)
model.seats.pvda.int.long <- model.seats.pvda.int %>%
  pivot_longer(cols = starts_with("fit"), names_to = "level", values_to = "probability")

# Define labels for the party categories used in the plot
party.labels <- c("Other", "PVDA")

# Convert 'level' column into a factor with a specified order for plotting
model.seats.pvda.int.long$level <- factor(model.seats.pvda.int.long$level, 
                                          levels = c("fit.Fewer", "fit.Same", "fit.More"))

# Open a PNG device to save the plot to file with specified resolution and dimensions
png("model_seats_pvda_int_long_g.png", units="in", width=10, height=7, res=1200)

# Create the plot
# --- Plotting: figure creation starts here. Review aesthetics (aes) and any facet/scales for readability.
ggplot(model.seats.pvda.int.long, aes(x = type, y = probability, color = party_pvda, linetype = party_pvda)) +
  geom_point() +  # Add points for each prediction
  # Add lines for each outcome level (fewer, same, more)
  geom_line(data = subset(model.seats.pvda.int.long, level == "fit.Fewer"),
            aes(x = type, y = probability, group = party_pvda, color = party_pvda, linetype = party_pvda)) +
  geom_line(data = subset(model.seats.pvda.int.long, level == "fit.Same"),
            aes(x = type, y = probability, group = party_pvda, color = party_pvda, linetype = party_pvda)) +
  geom_line(data = subset(model.seats.pvda.int.long, level == "fit.More"),
            aes(x = type, y = probability, group = party_pvda, color = party_pvda, linetype = party_pvda)) +
  # Split the plot into separate facets for each outcome category
  facet_grid(survey ~ level, 
             labeller = labeller(
               level = c("fit.Fewer" = "Fewer seats*", 
                         "fit.Same" = "Same as now*", 
                         "fit.More" = "More seats*"))) +
  # Axis and legend labels
  labs(title = "",
       x = "",
       y = "Predicted Probability",
       color = "Party Preference") +
  # Y-axis: fixed scale from 0 to 1
  scale_y_continuous(limits = c(0,1), breaks = seq(0, 1, 0.2)) +
  # X-axis: rename values 0 and 1 to descriptive labels
  scale_x_discrete(labels = c("0" = "Citizens", 
                              "1" = "Politicians")) +
  # Apply a minimal theme and customize appearance
  theme_minimal() +
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = 1),
    panel.spacing.x = unit(1.5, "lines"),  # Space between facets
    panel.spacing.y = unit(1.5, "lines"),  # Space between facets
    axis.title.x = element_text(size = 17, margin = margin(t = 10)),  # X-axis title styling
    axis.title.y = element_text(size = 17, margin = margin(r = 12)),  # Y-axis title styling
    axis.text.x = element_text(size = 15),    # X-axis tick labels
    axis.text.y = element_text(size = 15),    # Y-axis tick labels
    strip.text = element_text(size = 16),     # Facet label styling
    legend.position = "bottom",               # Move legend to bottom
    legend.title = element_text(size = 15),
    legend.text = element_text(size = 14),
    legend.key.height = unit(.8, "lines"),
    legend.key.width = unit(1, 'cm')) +
  # Format the legend: 2 columns, manual color and line type mappings
  guides(color = guide_legend(ncol = 2)) +
  scale_color_manual("Party Preference", values = c("#636363", "#990000"), labels = party.labels.pvda) +
  scale_linetype_manual("Party Preference", values = c("dashed", "solid"), labels = party.labels.pvda)

# Close the graphics device to save the plot
dev.off()

# =============================================================================
# 33. Predicted probabilities: Getting more or less seats (Groen)
# =============================================================================

# Predict class probabilities using the trained model for Groen
model.seats.groen.pr <- predict(model.seats_groen, newdata = model.groen.df, 
                                type = "prob")

# Combine the original data with predicted probabilities
model.seats.groen.int <- cbind(model.groen.df, model.seats.groen.pr)

# Reshape data from wide to long format for ggplot
model.seats.groen.int.long <- model.seats.groen.int %>%
  pivot_longer(cols = starts_with("fit"), names_to = "level", values_to = "probability")

# Define labels for the party categories used in the plot
party.labels.groen <- c("Other", "Groen")

# Convert 'level' column into a factor with a specified order for plotting
model.seats.groen.int.long$level <- factor(model.seats.groen.int.long$level, 
                                           levels = c("fit.Fewer", "fit.Same", "fit.More"))

# Save the plot
png("model_seats_groen_int_long_g.png", units="in", width=10, height=7, res=1200)

# --- Plotting: figure creation starts here. Review aesthetics (aes) and any facet/scales for readability.
ggplot(model.seats.groen.int.long, aes(x = type, y = probability, color = party_groen, linetype = party_groen)) +
  geom_point() +
  geom_line(data = subset(model.seats.groen.int.long, level == "fit.Fewer"),
            aes(group = party_groen)) +
  geom_line(data = subset(model.seats.groen.int.long, level == "fit.Same"),
            aes(group = party_groen)) +
  geom_line(data = subset(model.seats.groen.int.long, level == "fit.More"),
            aes(group = party_groen)) +
  facet_grid(survey ~ level, 
             labeller = labeller(
               level = c("fit.Fewer" = "Fewer seats", 
                         "fit.Same" = "Same as now", 
                         "fit.More" = "More seats"))) +
  labs(title = "",
       x = "",
       y = "Predicted Probability",
       color = "Party Preference") +
  scale_y_continuous(limits = c(0,1), breaks = seq(0, 1, 0.2)) +
  scale_x_discrete(labels = c("0" = "Citizens", "1" = "Politicians")) +
  theme_minimal() +
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = 1),
    panel.spacing.x = unit(1.5, "lines"),  # Space between facets
    panel.spacing.y = unit(1.5, "lines"),  # Space between facets
    axis.title.x = element_text(size = 17, margin = margin(t = 10)),
    axis.title.y = element_text(size = 17, margin = margin(r = 12)),
    axis.text.x = element_text(size = 15),
    axis.text.y = element_text(size = 15),
    strip.text = element_text(size = 16),
    legend.position = "bottom",
    legend.title = element_text(size = 15),
    legend.text = element_text(size = 14),
    legend.key.height = unit(.8, "lines"),
    legend.key.width = unit(1, 'cm')) +
  guides(color = guide_legend(ncol = 2)) +
  scale_color_manual("Party Preference", values = c("#636363", "#04847A"), labels = party.labels.groen) +
  scale_linetype_manual("Party Preference", values = c("dashed", "solid"), labels = party.labels.groen)

dev.off()

# =============================================================================
# 34. Predicted probabilities: Getting more or less seats (Vooruit)
# =============================================================================

# Predict class probabilities using the trained model for Vooruit
model.seats.vooruit.pr <- predict(model.seats_vooruit, newdata = model.vooruit.df, 
                                  type = "prob")  # Return probabilities for each class

# Combine the original data with predicted probabilities
model.seats.vooruit.int <- cbind(model.vooruit.df, model.seats.vooruit.pr)

# Reshape data from wide to long format for ggplot
model.seats.vooruit.int.long <- model.seats.vooruit.int %>%
  pivot_longer(cols = starts_with("fit"), names_to = "level", values_to = "probability")

# Define labels for the party categories used in the plot
party.labels.vooruit <- c("Other", "Vooruit")

# Convert 'level' column into a factor with a specified order
model.seats.vooruit.int.long$level <- factor(model.seats.vooruit.int.long$level, 
                                             levels = c("fit.Fewer", "fit.Same", "fit.More"))

# Open a PNG device to save the plot to file
png("model_seats_vooruit_int_long_g.png", units="in", width=10, height=7, res=1200)

# Create the plot
# --- Plotting: figure creation starts here. Review aesthetics (aes) and any facet/scales for readability.
ggplot(model.seats.vooruit.int.long, aes(x = type, y = probability, color = party_vooruit, linetype = party_vooruit)) +
  geom_point() +
  geom_line(data = subset(model.seats.vooruit.int.long, level == "fit.Fewer"),
            aes(group = party_vooruit)) +
  geom_line(data = subset(model.seats.vooruit.int.long, level == "fit.Same"),
            aes(group = party_vooruit)) +
  geom_line(data = subset(model.seats.vooruit.int.long, level == "fit.More"),
            aes(group = party_vooruit)) +
  facet_grid(survey ~ level, 
             labeller = labeller(
               level = c("fit.Fewer" = "Fewer seats*", 
                         "fit.Same" = "Same as now", 
                         "fit.More" = "More seats"))) +
  labs(title = "",
       x = "",
       y = "Predicted Probability",
       color = "Party Preference") +
  scale_y_continuous(limits = c(0,1), breaks = seq(0, 1, 0.2)) +
  scale_x_discrete(labels = c("0" = "Citizens", "1" = "Politicians")) +
  theme_minimal() +
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = 1),
    panel.spacing.x = unit(1.5, "lines"),  # Space between facets
    panel.spacing.y = unit(1.5, "lines"),  # Space between facets
    axis.title.x = element_text(size = 17, margin = margin(t = 10)),
    axis.title.y = element_text(size = 17, margin = margin(r = 12)),
    axis.text.x = element_text(size = 15),
    axis.text.y = element_text(size = 15),
    strip.text = element_text(size = 16),
    legend.position = "bottom",
    legend.title = element_text(size = 15),
    legend.text = element_text(size = 14),
    legend.key.height = unit(.8, "lines"),
    legend.key.width = unit(1, 'cm')) +
  guides(color = guide_legend(ncol = 2)) +
  scale_color_manual("Party Preference", values = c("#636363", "#FF0000"), labels = party.labels.vooruit) +
  scale_linetype_manual("Party Preference", values = c("dashed", "solid"), labels = party.labels.vooruit)

# Close the graphics device to save the plot
dev.off()

# =============================================================================
# 35. Predicted probabilities: Getting more or less seats (CD&V)
# =============================================================================

model.seats.cdv.pr <- predict(model.seats_cdv, newdata = model.cdv.df, type = "prob")
model.seats.cdv.int <- cbind(model.cdv.df, model.seats.cdv.pr)

model.seats.cdv.int.long <- model.seats.cdv.int %>%
  pivot_longer(cols = starts_with("fit"), names_to = "level", values_to = "probability")

party.labels.cdv <- c("Other", "CD&V")

model.seats.cdv.int.long$level <- factor(model.seats.cdv.int.long$level, 
                                         levels = c("fit.Fewer", "fit.Same", "fit.More"))

png("model_seats_cdv_int_long_g.png", units="in", width=10, height=7, res=1200)

# --- Plotting: figure creation starts here. Review aesthetics (aes) and any facet/scales for readability.
ggplot(model.seats.cdv.int.long, aes(x = type, y = probability, color = party_cdv, linetype = party_cdv)) +
  geom_point() +
  geom_line(data = subset(model.seats.cdv.int.long, level == "fit.Fewer"), aes(group = party_cdv)) +
  geom_line(data = subset(model.seats.cdv.int.long, level == "fit.Same"), aes(group = party_cdv)) +
  geom_line(data = subset(model.seats.cdv.int.long, level == "fit.More"), aes(group = party_cdv)) +
  facet_grid(survey ~ level, 
             labeller = labeller(
               level = c("fit.Fewer" = "Fewer seats*", 
                         "fit.Same" = "Same as now*", 
                         "fit.More" = "More seats*"))) +
  labs(title = "",
       x = "",
       y = "Predicted Probability",
       color = "Party Preference") +
  scale_y_continuous(limits = c(0,1), breaks = seq(0, 1, 0.2)) +
  scale_x_discrete(labels = c("0" = "Citizens", "1" = "Politicians")) +
  theme_minimal() +
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = 1),
    panel.spacing.x = unit(1.5, "lines"),  # Space between facets
    panel.spacing.y = unit(1.5, "lines"),  # Space between facets
    axis.title.x = element_text(size = 17, margin = margin(t = 10)),
    axis.title.y = element_text(size = 17, margin = margin(r = 12)),
    axis.text.x = element_text(size = 15),
    axis.text.y = element_text(size = 15),
    strip.text = element_text(size = 16),
    legend.position = "bottom",
    legend.title = element_text(size = 15),
    legend.text = element_text(size = 14),
    legend.key.height = unit(.8, "lines"),
    legend.key.width = unit(1, 'cm')) +
  guides(color = guide_legend(ncol = 2)) +
  scale_color_manual("Party Preference", values = c("#636363", "#FF6000"), labels = party.labels.cdv) +
  scale_linetype_manual("Party Preference", values = c("dashed", "solid"), labels = party.labels.cdv)

dev.off()

# =============================================================================
# 36. Predicted probabilities: Getting more or less seats (Open Vld)
# =============================================================================

model.seats.ovld.pr <- predict(model.seats_ovld, newdata = model.ovld.df, type = "prob")
model.seats.ovld.int <- cbind(model.ovld.df, model.seats.ovld.pr)

model.seats.ovld.int.long <- model.seats.ovld.int %>%
  pivot_longer(cols = starts_with("fit"), names_to = "level", values_to = "probability")

party.labels.ovld <- c("Other", "Open Vld")

model.seats.ovld.int.long$level <- factor(model.seats.ovld.int.long$level, 
                                          levels = c("fit.Fewer", "fit.Same", "fit.More"))

png("model_seats_ovld_int_long_g.png", units="in", width=10, height=7, res=1200)

# --- Plotting: figure creation starts here. Review aesthetics (aes) and any facet/scales for readability.
ggplot(model.seats.ovld.int.long, aes(x = type, y = probability, color = party_ovld, linetype = party_ovld)) +
  geom_point() +
  geom_line(data = subset(model.seats.ovld.int.long, level == "fit.Fewer"), aes(group = party_ovld)) +
  geom_line(data = subset(model.seats.ovld.int.long, level == "fit.Same"), aes(group = party_ovld)) +
  geom_line(data = subset(model.seats.ovld.int.long, level == "fit.More"), aes(group = party_ovld)) +
  facet_grid(survey ~ level, 
             labeller = labeller(
               level = c("fit.Fewer" = "Fewer seats*", 
                         "fit.Same" = "Same as now*", 
                         "fit.More" = "More seats*"))) +
  labs(title = "",
       x = "",
       y = "Predicted Probability",
       color = "Party Preference") +
  scale_y_continuous(limits = c(0,1), breaks = seq(0, 1, 0.2)) +
  scale_x_discrete(labels = c("0" = "Citizens", "1" = "Politicians")) +
  theme_minimal() +
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = 1),
    panel.spacing.x = unit(1.5, "lines"),  # Space between facets
    panel.spacing.y = unit(1.5, "lines"),  # Space between facets
    axis.title.x = element_text(size = 17, margin = margin(t = 10)),
    axis.title.y = element_text(size = 17, margin = margin(r = 12)),
    axis.text.x = element_text(size = 15),
    axis.text.y = element_text(size = 15),
    strip.text = element_text(size = 16),
    legend.position = "bottom",
    legend.title = element_text(size = 15),
    legend.text = element_text(size = 14),
    legend.key.height = unit(.8, "lines"),
    legend.key.width = unit(1, 'cm')) +
  guides(color = guide_legend(ncol = 2)) +
  scale_color_manual("Party Preference", values = c("#636363", "#005DAA"), labels = party.labels.ovld) +
  scale_linetype_manual("Party Preference", values = c("dashed", "solid"), labels = party.labels.ovld)

dev.off()

# =============================================================================
# 37. Predicted probabilities: Getting more or less seats (N-VA)
# =============================================================================

model.seats.nva.pr <- predict(model.seats_nva, newdata = model.nva.df, type = "prob")
model.seats.nva.int <- cbind(model.nva.df, model.seats.nva.pr)

model.seats.nva.int.long <- model.seats.nva.int %>%
  pivot_longer(cols = starts_with("fit"), names_to = "level", values_to = "probability")

party.labels.nva <- c("Other", "N-VA")

model.seats.nva.int.long$level <- factor(model.seats.nva.int.long$level, 
                                         levels = c("fit.Fewer", "fit.Same", "fit.More"))

png("model_seats_nva_int_long_g.png", units="in", width=10, height=7, res=1200)

# --- Plotting: figure creation starts here. Review aesthetics (aes) and any facet/scales for readability.
ggplot(model.seats.nva.int.long, aes(x = type, y = probability, color = party_nva, linetype = party_nva)) +
  geom_point() +
  geom_line(data = subset(model.seats.nva.int.long, level == "fit.Fewer"), aes(group = party_nva)) +
  geom_line(data = subset(model.seats.nva.int.long, level == "fit.Same"), aes(group = party_nva)) +
  geom_line(data = subset(model.seats.nva.int.long, level == "fit.More"), aes(group = party_nva)) +
  facet_grid(survey ~ level, 
             labeller = labeller(
               level = c("fit.Fewer" = "Fewer seats", 
                         "fit.Same" = "Same as now*", 
                         "fit.More" = "More seats"))) +
  labs(title = "",
       x = "",
       y = "Predicted Probability",
       color = "Party Preference") +
  scale_y_continuous(limits = c(0,1), breaks = seq(0, 1, 0.2)) +
  scale_x_discrete(labels = c("0" = "Citizens", "1" = "Politicians")) +
  theme_minimal() +
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = 1),
    panel.spacing.x = unit(1.5, "lines"),  # Space between facets
    panel.spacing.y = unit(1.5, "lines"),  # Space between facets
    axis.title.x = element_text(size = 17, margin = margin(t = 10)),
    axis.title.y = element_text(size = 17, margin = margin(r = 12)),
    axis.text.x = element_text(size = 15),
    axis.text.y = element_text(size = 15),
    strip.text = element_text(size = 16),
    legend.position = "bottom",
    legend.title = element_text(size = 15),
    legend.text = element_text(size = 14),
    legend.key.height = unit(.8, "lines"),
    legend.key.width = unit(1, 'cm')) +
  guides(color = guide_legend(ncol = 2)) +
  scale_color_manual("Party Preference", values = c("#636363", "#FCBD1B"), labels = party.labels.nva) +
  scale_linetype_manual("Party Preference", values = c("dashed", "solid"), labels = party.labels.nva)

dev.off()

# =============================================================================
# 38. Predicted probabilities: Getting more or less seats (Vlaams Belang)
# =============================================================================

model.seats.vb.pr <- predict(model.seats_vb, newdata = model.vb.df, type = "prob")
model.seats.vb.int <- cbind(model.vb.df, model.seats.vb.pr)

model.seats.vb.int.long <- model.seats.vb.int %>%
  pivot_longer(cols = starts_with("fit"), names_to = "level", values_to = "probability")

party.labels.vb <- c("Other", "VB")

model.seats.vb.int.long$level <- factor(model.seats.vb.int.long$level, 
                                        levels = c("fit.Fewer", "fit.Same", "fit.More"))

png("model_seats_vb_int_long_g.png", units="in", width=10, height=7, res=1200)

# --- Plotting: figure creation starts here. Review aesthetics (aes) and any facet/scales for readability.
ggplot(model.seats.vb.int.long, aes(x = type, y = probability, color = party_vb, linetype = party_vb)) +
  geom_point() +
  geom_line(data = subset(model.seats.vb.int.long, level == "fit.Fewer"), aes(group = party_vb)) +
  geom_line(data = subset(model.seats.vb.int.long, level == "fit.Same"), aes(group = party_vb)) +
  geom_line(data = subset(model.seats.vb.int.long, level == "fit.More"), aes(group = party_vb)) +
  facet_grid(survey ~ level, 
             labeller = labeller(
               level = c("fit.Fewer" = "Fewer seats*", 
                         "fit.Same" = "Same as now*", 
                         "fit.More" = "More seats*"))) +
  labs(title = "",
       x = "",
       y = "Predicted Probability",
       color = "Party Preference") +
  scale_y_continuous(limits = c(0,1), breaks = seq(0, 1, 0.2)) +
  scale_x_discrete(labels = c("0" = "Citizens", "1" = "Politicians")) +
  theme_minimal() +
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = 1),
    panel.spacing.x = unit(1.5, "lines"),  # Space between facets
    panel.spacing.y = unit(1.5, "lines"),  # Space between facets
    axis.title.x = element_text(size = 17, margin = margin(t = 10)),
    axis.title.y = element_text(size = 17, margin = margin(r = 12)),
    axis.text.x = element_text(size = 15),
    axis.text.y = element_text(size = 15),
    strip.text = element_text(size = 16),
    legend.position = "bottom",
    legend.title = element_text(size = 15),
    legend.text = element_text(size = 14),
    legend.key.height = unit(.8, "lines"),
    legend.key.width = unit(1, 'cm')) +
  guides(color = guide_legend(ncol = 2)) +
  scale_color_manual("Party Preference", values = c("#636363", "black"), labels = party.labels.vb) +
  scale_linetype_manual("Party Preference", values = c("dashed", "solid"), labels = party.labels.vb)

dev.off()

# =============================================================================
# 39. Structure predicted probabilities into wide format (PVDA)
# =============================================================================

# Copy the long-format predictions into a new object for PVDA
prob.pvda <- model.seats.pvda.int.long

# Select relevant columns: party, respondent type (citizen/politician), outcome level, and predicted probability
# --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
prob.pvda <- prob.pvda %>% dplyr::select(party_pvda, type, survey, level, probability)

# Rename 'party_pvda' to generic 'party' to standardize column names across parties
prob.pvda <- prob.pvda %>% dplyr::rename(party = party_pvda)

# Sort the data by outcome level (fit.Fewer, fit.Same, fit.More) for clarity
prob.pvda <- prob.pvda[order(prob.pvda$level),]

# Separate data by respondent type:
# cf = citizens (type == 0), pf = politicians (type == 1)
# --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
cf.prob.pvda <- prob.pvda %>% filter(type == 0)
# --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
pf.prob.pvda <- prob.pvda %>% filter(type == 1)

# Reshape citizen data from long to wide:
# Each outcome level becomes a separate column (fit.Fewer, fit.Same, fit.More)
cf.prob.pvda <- cf.prob.pvda %>%
  pivot_wider(names_from = level, values_from = probability)

# Further reshape so each value is nested under its corresponding party label
cf.prob.pvda <- cf.prob.pvda %>%
  pivot_wider(names_from = party, values_from = c(fit.Fewer, fit.Same, fit.More))

# Add a 'party' column to identify this row as belonging to PVDA
cf.prob.pvda$party <- "PVDA"

# Repeat the same reshaping process for politician data
pf.prob.pvda <- pf.prob.pvda %>%
  pivot_wider(names_from = level, values_from = probability)

pf.prob.pvda <- pf.prob.pvda %>%
  pivot_wider(names_from = party, values_from = c(fit.Fewer, fit.Same, fit.More))

pf.prob.pvda$party <- "PVDA"

# Merge reshaped citizen and politician data by party
# This aligns each party's citizen and politician predictions side by side
prob.pvda <- merge(cf.prob.pvda, pf.prob.pvda, by = c("party", "survey"), all = FALSE)

# =============================================================================
# 40. Structure predicted probabilities into wide format (Groen)
# =============================================================================

# Copy the long-format predictions into a new object for Groen
prob.groen <- model.seats.groen.int.long

# Select relevant columns: party, respondent type (citizen/politician), outcome level, and predicted probability
# --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
prob.groen <- prob.groen %>% dplyr::select(party_groen, type, survey, level, probability)

# Rename 'party_groen' to generic 'party' to standardize column names across parties
prob.groen <- prob.groen %>% dplyr::rename(party = party_groen)

# Sort the data by outcome level (fit.Fewer, fit.Same, fit.More) for clarity
prob.groen <- prob.groen[order(prob.groen$level),]

# Separate data by respondent type:
# cf = citizens (type == 0), pf = politicians (type == 1)
# --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
cf.prob.groen <- prob.groen %>% filter(type == 0)
# --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
pf.prob.groen <- prob.groen %>% filter(type == 1)

# Reshape citizen data from long to wide:
# Each outcome level becomes a separate column (fit.Fewer, fit.Same, fit.More)
cf.prob.groen <- cf.prob.groen %>%
  pivot_wider(names_from = level, values_from = probability)

# Further reshape so each value is nested under its corresponding party label
cf.prob.groen <- cf.prob.groen %>%
  pivot_wider(names_from = party, values_from = c(fit.Fewer, fit.Same, fit.More))

# Add a 'party' column to identify this row as belonging to Groen
cf.prob.groen$party <- "Groen"

# Repeat the same reshaping process for politician data
pf.prob.groen <- pf.prob.groen %>%
  pivot_wider(names_from = level, values_from = probability)

pf.prob.groen <- pf.prob.groen %>%
  pivot_wider(names_from = party, values_from = c(fit.Fewer, fit.Same, fit.More))

pf.prob.groen$party <- "Groen"

# Merge reshaped citizen and politician data by party
# This aligns each party's citizen and politician predictions side by side
prob.groen <- merge(cf.prob.groen, pf.prob.groen, by = c("party", "survey"), all = FALSE)

# =============================================================================
# 41. Structure predicted probabilities into wide format (Vooruit)
# =============================================================================

# Copy the long-format predictions into a new object for Vooruit
prob.vooruit <- model.seats.vooruit.int.long

# Select relevant columns: party, respondent type (citizen/politician), outcome level, and predicted probability
# --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
prob.vooruit <- prob.vooruit %>% dplyr::select(party_vooruit, type, survey, level, probability)

# Rename 'party_vooruit' to generic 'party' to standardize column names across parties
prob.vooruit <- prob.vooruit %>% dplyr::rename(party = party_vooruit)

# Sort the data by outcome level (fit.Fewer, fit.Same, fit.More) for clarity
prob.vooruit <- prob.vooruit[order(prob.vooruit$level),]

# Separate data by respondent type:
# cf = citizens (type == 0), pf = politicians (type == 1)
# --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
cf.prob.vooruit <- prob.vooruit %>% filter(type == 0)
# --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
pf.prob.vooruit <- prob.vooruit %>% filter(type == 1)

# Reshape citizen data from long to wide:
# Each outcome level becomes a separate column (fit.Fewer, fit.Same, fit.More)
cf.prob.vooruit <- cf.prob.vooruit %>%
  pivot_wider(names_from = level, values_from = probability)

# Further reshape so each value is nested under its corresponding party label
cf.prob.vooruit <- cf.prob.vooruit %>%
  pivot_wider(names_from = party, values_from = c(fit.Fewer, fit.Same, fit.More))

# Add a 'party' column to identify this row as belonging to Vooruit
cf.prob.vooruit$party <- "Vooruit"

# Repeat the same reshaping process for politician data
pf.prob.vooruit <- pf.prob.vooruit %>%
  pivot_wider(names_from = level, values_from = probability)

pf.prob.vooruit <- pf.prob.vooruit %>%
  pivot_wider(names_from = party, values_from = c(fit.Fewer, fit.Same, fit.More))

pf.prob.vooruit$party <- "Vooruit"

# Merge reshaped citizen and politician data by party
# This aligns each party's citizen and politician predictions side by side
prob.vooruit <- merge(cf.prob.vooruit, pf.prob.vooruit, by = c("party", "survey"), all = FALSE)

# =============================================================================
# 42. Structure predicted probabilities into wide format (CD&V)
# =============================================================================

# Copy the long-format predictions into a new object for CD&V
prob.cdv <- model.seats.cdv.int.long

# Select relevant columns: party, respondent type (citizen/politician), outcome level, and predicted probability
# --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
prob.cdv <- prob.cdv %>% dplyr::select(party_cdv, type, survey, level, probability)

# Rename 'party_cdv' to generic 'party' to standardize column names across parties
prob.cdv <- prob.cdv %>% dplyr::rename(party = party_cdv)

# Sort the data by outcome level (fit.Fewer, fit.Same, fit.More) for clarity
prob.cdv <- prob.cdv[order(prob.cdv$level),]

# Separate data by respondent type:
# cf = citizens (type == 0), pf = politicians (type == 1)
# --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
cf.prob.cdv <- prob.cdv %>% filter(type == 0)
# --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
pf.prob.cdv <- prob.cdv %>% filter(type == 1)

# Reshape citizen data from long to wide:
# Each outcome level becomes a separate column (fit.Fewer, fit.Same, fit.More)
cf.prob.cdv <- cf.prob.cdv %>%
  pivot_wider(names_from = level, values_from = probability)

# Further reshape so each value is nested under its corresponding party label
cf.prob.cdv <- cf.prob.cdv %>%
  pivot_wider(names_from = party, values_from = c(fit.Fewer, fit.Same, fit.More))

# Add a 'party' column to identify this row as belonging to CD&V
cf.prob.cdv$party <- "CD&V"

# Repeat the same reshaping process for politician data
pf.prob.cdv <- pf.prob.cdv %>%
  pivot_wider(names_from = level, values_from = probability)

pf.prob.cdv <- pf.prob.cdv %>%
  pivot_wider(names_from = party, values_from = c(fit.Fewer, fit.Same, fit.More))

pf.prob.cdv$party <- "CD&V"

# Merge reshaped citizen and politician data by party
# This aligns each party's citizen and politician predictions side by side
prob.cdv <- merge(cf.prob.cdv, pf.prob.cdv, by = c("party", "survey"), all = FALSE)

# =============================================================================
# 43. Structure predicted probabilities into wide format (Open Vld)
# =============================================================================

# Copy the long-format predictions into a new object for Open Vld
prob.ovld <- model.seats.ovld.int.long

# Select relevant columns: party, respondent type (citizen/politician), outcome level, and predicted probability
# --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
prob.ovld <- prob.ovld %>% dplyr::select(party_ovld, type, survey, level, probability)

# Rename 'party_ovld' to generic 'party' to standardize column names across parties
prob.ovld <- prob.ovld %>% dplyr::rename(party = party_ovld)

# Sort the data by outcome level (fit.Fewer, fit.Same, fit.More) for clarity
prob.ovld <- prob.ovld[order(prob.ovld$level),]

# Separate data by respondent type:
# cf = citizens (type == 0), pf = politicians (type == 1)
# --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
cf.prob.ovld <- prob.ovld %>% filter(type == 0)
# --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
pf.prob.ovld <- prob.ovld %>% filter(type == 1)

# Reshape citizen data from long to wide:
# Each outcome level becomes a separate column (fit.Fewer, fit.Same, fit.More)
cf.prob.ovld <- cf.prob.ovld %>%
  pivot_wider(names_from = level, values_from = probability)

# Further reshape so each value is nested under its corresponding party label
cf.prob.ovld <- cf.prob.ovld %>%
  pivot_wider(names_from = party, values_from = c(fit.Fewer, fit.Same, fit.More))

# Add a 'party' column to identify this row as belonging to Open Vld
cf.prob.ovld$party <- "Open Vld"

# Repeat the same reshaping process for politician data
pf.prob.ovld <- pf.prob.ovld %>%
  pivot_wider(names_from = level, values_from = probability)

pf.prob.ovld <- pf.prob.ovld %>%
  pivot_wider(names_from = party, values_from = c(fit.Fewer, fit.Same, fit.More))

pf.prob.ovld$party <- "Open Vld"

# Merge reshaped citizen and politician data by party
# This aligns each party's citizen and politician predictions side by side
prob.ovld <- merge(cf.prob.ovld, pf.prob.ovld, by = c("party", "survey"), all = FALSE)

# =============================================================================
# 44. Structure predicted probabilities into wide format (N-VA)
# =============================================================================

# Copy the long-format predictions into a new object for N-VA
prob.nva <- model.seats.nva.int.long

# Select relevant columns: party, respondent type (citizen/politician), outcome level, and predicted probability
# --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
prob.nva <- prob.nva %>% dplyr::select(party_nva, type, survey, level, probability)

# Rename 'party_nva' to generic 'party' to standardize column names across parties
prob.nva <- prob.nva %>% dplyr::rename(party = party_nva)

# Sort the data by outcome level (fit.Fewer, fit.Same, fit.More) for clarity
prob.nva <- prob.nva[order(prob.nva$level),]

# Separate data by respondent type:
# cf = citizens (type == 0), pf = politicians (type == 1)
# --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
cf.prob.nva <- prob.nva %>% filter(type == 0)
# --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
pf.prob.nva <- prob.nva %>% filter(type == 1)

# Reshape citizen data from long to wide:
# Each outcome level becomes a separate column (fit.Fewer, fit.Same, fit.More)
cf.prob.nva <- cf.prob.nva %>%
  pivot_wider(names_from = level, values_from = probability)

# Further reshape so each value is nested under its corresponding party label
cf.prob.nva <- cf.prob.nva %>%
  pivot_wider(names_from = party, values_from = c(fit.Fewer, fit.Same, fit.More))

# Add a 'party' column to identify this row as belonging to N-VA
cf.prob.nva$party <- "N-VA"

# Repeat the same reshaping process for politician data
pf.prob.nva <- pf.prob.nva %>%
  pivot_wider(names_from = level, values_from = probability)

pf.prob.nva <- pf.prob.nva %>%
  pivot_wider(names_from = party, values_from = c(fit.Fewer, fit.Same, fit.More))

pf.prob.nva$party <- "N-VA"

# Merge reshaped citizen and politician data by party
# This aligns each party's citizen and politician predictions side by side
prob.nva <- merge(cf.prob.nva, pf.prob.nva, by = c("party", "survey"), all = FALSE)

# =============================================================================
# 45. Structure predicted probabilities into wide format (Vlaams Belang)
# =============================================================================

# Copy the long-format predictions into a new object for VLAAMS BELANG
prob.vb <- model.seats.vb.int.long

# Select relevant columns: party, respondent type (citizen/politician), outcome level, and predicted probability
# --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
prob.vb <- prob.vb %>% dplyr::select(party_vb, type, survey, level, probability)

# Rename 'party_vb' to generic 'party' to standardize column names across parties
prob.vb <- prob.vb %>% dplyr::rename(party = party_vb)

# Sort the data by outcome level (fit.Fewer, fit.Same, fit.More) for clarity
prob.vb <- prob.vb[order(prob.vb$level),]

# Separate data by respondent type:
# cf = citizens (type == 0), pf = politicians (type == 1)
# --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
cf.prob.vb <- prob.vb %>% filter(type == 0)
# --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
pf.prob.vb <- prob.vb %>% filter(type == 1)

# Reshape citizen data from long to wide:
# Each outcome level becomes a separate column (fit.Fewer, fit.Same, fit.More)
cf.prob.vb <- cf.prob.vb %>%
  pivot_wider(names_from = level, values_from = probability)

# Further reshape so each value is nested under its corresponding party label
cf.prob.vb <- cf.prob.vb %>%
  pivot_wider(names_from = party, values_from = c(fit.Fewer, fit.Same, fit.More))

# Add a 'party' column to identify this row as belonging to VLAAMS BELANG
cf.prob.vb$party <- "Vlaams Belang"

# Repeat the same reshaping process for politician data
pf.prob.vb <- pf.prob.vb %>%
  pivot_wider(names_from = level, values_from = probability)

pf.prob.vb <- pf.prob.vb %>%
  pivot_wider(names_from = party, values_from = c(fit.Fewer, fit.Same, fit.More))

pf.prob.vb$party <- "Vlaams Belang"

# Merge reshaped citizen and politician data by party
# This aligns each party's citizen and politician predictions side by side
prob.vb <- merge(cf.prob.vb, pf.prob.vb, by = c("party", "survey"), all = FALSE)

################################################################################
## Combine, calculate, and visualize predicted probability differences
## for multiple parties (citizens vs elites)
################################################################################

# Combine probability data frames for all parties into one
prob <- rbind(prob.pvda, prob.groen, prob.vooruit, prob.cdv, prob.ovld, prob.nva, prob.vb)

# Calculate the absolute difference in predicted probabilities (scaled by 100) between citizens (type 0)
# and politicians (type 1) for 'fewer', 'same', and 'more' seat changes - mass (citizen) perspective
prob$mass.fewer <- abs(prob$fit.Fewer_0.x - prob$fit.Fewer_1.x) * 100
prob$mass.same <- abs(prob$fit.Same_0.x - prob$fit.Same_1.x) * 100
prob$mass.more <- abs(prob$fit.More_0.x - prob$fit.More_1.x) * 100

# Similarly, calculate the absolute differences for the elite (politician) perspective
prob$elite.fewer <- abs(prob$fit.Fewer_0.y - prob$fit.Fewer_1.y) * 100
prob$elite.same <- abs(prob$fit.Same_0.y - prob$fit.Same_1.y) * 100
prob$elite.more <- abs(prob$fit.More_0.y - prob$fit.More_1.y) * 100

# Compute the absolute differences between mass and elite differences for each outcome
prob$diff.fewer <- abs(prob$mass.fewer - prob$elite.fewer)
prob$diff.same <- abs(prob$mass.same - prob$elite.same)
prob$diff.more <- abs(prob$mass.more - prob$elite.more)

# Select and reorder columns to keep relevant variables for output
# --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
prob <- prob %>% dplyr::select(party, survey, mass.fewer, elite.fewer, diff.fewer,
                               mass.same, elite.same, diff.same,
                               mass.more, elite.more, diff.more)

# Export the summarized differences to an Excel file for further review
write_xlsx(prob, "prob.xlsx")

# =============================================================================
# 46. Prepare elite (politician) data for long-form plotting
# =============================================================================

# Select elite probabilities and party identifiers for long-format transformation
# --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
prob.elite <- prob %>% dplyr::select(survey, matches("elite|party", ignore.case = TRUE))

# Convert wide elite data to long format: one row per party and change type (fewer/same/more)
prob.elite.long <- prob.elite %>%
  pivot_longer(cols = starts_with("elite"), 
               names_to = "change", 
               names_prefix = "(.+?)\\.",  # Remove 'fewer.', 'same.', 'more.' prefix from names
               values_to = "value")

# Add a type indicator "P" for politicians (elite)
prob.elite.long$type <- "P"

# Similarly, prepare the elite difference data for plotting
# --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
prob.elite.diff <- prob %>% dplyr::select(survey, matches("diff|party", ignore.case = TRUE))

prob.elite.diff.long <- prob.elite.diff %>%
  pivot_longer(cols = starts_with("diff"), 
               names_to = "change", 
               names_prefix = "(.+?)\\.", 
               values_to = "diff")

prob.elite.diff.long$type <- "P"

# Merge elite values and their corresponding differences by party, change, and type
prob.elite.df <- merge(prob.elite.long, prob.elite.diff.long, by = c("party", "change", "type", "survey"))

# Assign significance labels based on party and change type (pre-determined results)
prob.elite.df <- prob.elite.df %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  mutate(sig = case_when(
    party == "PVDA" & change == "fewer" & survey == "2023" ~ "Significant",
    party == "PVDA" & change == "fewer" & survey == "2024" ~ "Significant",
    
    party == "Groen" & change == "fewer" & survey == "2023" ~ "Non-significant",
    party == "Groen" & change == "fewer" & survey == "2024" ~ "Non-significant",
    
    party == "Vooruit" & change == "fewer" & survey == "2023" ~ "Significant",
    party == "Vooruit" & change == "fewer" & survey == "2024" ~ "Significant",
    
    party == "CD&V" & change == "fewer" & survey == "2023" ~ "Non-significant",
    party == "CD&V" & change == "fewer" & survey == "2024" ~ "Non-significant",
    
    party == "Open Vld" & change == "fewer" & survey == "2023" ~ "Non-significant",
    party == "Open Vld" & change == "fewer" & survey == "2024" ~ "Non-significant",
    
    party == "N-VA" & change == "fewer" & survey == "2023" ~ "Non-significant",
    party == "N-VA" & change == "fewer" & survey == "2024" ~ "Non-significant",
    
    party == "Vlaams Belang" & change == "fewer" & survey == "2023" ~ "Significant",
    party == "Vlaams Belang" & change == "fewer" & survey == "2024" ~ "Significant",
    
    party == "PVDA" & change == "same" & survey == "2023" ~ "Significant",
    party == "PVDA" & change == "same" & survey == "2024" ~ "Significant",
    
    party == "Groen" & change == "same" & survey == "2023" ~ "Non-significant",
    party == "Groen" & change == "same" & survey == "2024" ~ "Non-significant",
    
    party == "Vooruit" & change == "same" & survey == "2023" ~ "Significant",
    party == "Vooruit" & change == "same" & survey == "2024" ~ "Significant",
    
    party == "CD&V" & change == "same" & survey == "2023" ~ "Significant",
    party == "CD&V" & change == "same" & survey == "2024" ~ "Non-significant",
    
    party == "Open Vld" & change == "same" & survey == "2023" ~ "Non-significant",
    party == "Open Vld" & change == "same" & survey == "2024" ~ "Non-significant",
    
    party == "N-VA" & change == "same" & survey == "2023" ~ "Non-significant",
    party == "N-VA" & change == "same" & survey == "2024" ~ "Non-significant",
    
    party == "Vlaams Belang" & change == "same" & survey == "2023" ~ "Non-significant",
    party == "Vlaams Belang" & change == "same" & survey == "2024" ~ "Non-significant",
    
    party == "PVDA" & change == "more" & survey == "2023" ~ "Significant",
    party == "PVDA" & change == "more" & survey == "2024" ~ "Significant",
    
    party == "Groen" & change == "more" & survey == "2023" ~ "Non-significant",
    party == "Groen" & change == "more" & survey == "2024" ~ "Non-significant",
    
    party == "Vooruit" & change == "more" & survey == "2023" ~ "Non-significant",
    party == "Vooruit" & change == "more" & survey == "2024" ~ "Non-significant",
    
    party == "CD&V" & change == "more" & survey == "2023" ~ "Non-significant",
    party == "CD&V" & change == "more" & survey == "2024" ~ "Non-significant",
    
    party == "Open Vld" & change == "more" & survey == "2023" ~ "Non-significant",
    party == "Open Vld" & change == "more" & survey == "2024" ~ "Non-significant",
    
    party == "N-VA" & change == "more" & survey == "2023" ~ "Significant",
    party == "N-VA" & change == "more" & survey == "2024" ~ "Non-significant",
    
    party == "Vlaams Belang" & change == "more" & survey == "2023" ~ "Significant",
    party == "Vlaams Belang" & change == "more" & survey == "2024" ~ "Significant"))

# Assign numeric IDs for parties to facilitate plotting order
prob.elite.df <- prob.elite.df %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  mutate(id = case_when(
    party == "PVDA" ~ 7,
    party == "Groen" ~ 6,
    party == "Vooruit" ~ 5,
    party == "CD&V" ~ 4,
    party == "Open Vld" ~ 3,
    party == "N-VA" ~ 2,
    party == "Vlaams Belang" ~ 1
  ))

# =============================================================================
# 47. Prepare mass (citizen) data for long-form plotting
# =============================================================================

# Select mass probabilities and party identifiers for long-format transformation
# --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
prob.mass <- prob %>% dplyr::select(survey, matches("mass|party", ignore.case = TRUE))

# Convert wide mass data to long format: one row per party and change type
prob.mass.long <- prob.mass %>%
  pivot_longer(cols = starts_with("mass"), 
               names_to = "change", 
               names_prefix = "(.+?)\\.",  # Strip prefix before change label
               values_to = "value")

# Add a type indicator "C" for citizens (mass)
prob.mass.long$type <- "C"

# Prepare difference values for mass data similarly
# --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
prob.mass.diff <- prob %>% dplyr::select(survey, matches("diff|party", ignore.case = TRUE))

prob.mass.diff.long <- prob.mass.diff %>%
  pivot_longer(cols = starts_with("diff"), 
               names_to = "change", 
               names_prefix = "(.+?)\\.", 
               values_to = "diff")

prob.mass.diff.long$type <- "C"

# Merge mass values and their corresponding differences by party, change, and type
prob.mass.df <- merge(prob.mass.long, prob.mass.diff.long, by = c("party", "change", "type", "survey"))

# Assign significance labels as before based on party and change type
prob.mass.df <- prob.mass.df %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  mutate(sig = case_when(
    party == "PVDA" & change == "fewer" & survey == "2023" ~ "Significant",
    party == "PVDA" & change == "fewer" & survey == "2024" ~ "Significant",
    
    party == "Groen" & change == "fewer" & survey == "2023" ~ "Non-significant",
    party == "Groen" & change == "fewer" & survey == "2024" ~ "Non-significant",
    
    party == "Vooruit" & change == "fewer" & survey == "2023" ~ "Significant",
    party == "Vooruit" & change == "fewer" & survey == "2024" ~ "Significant",
    
    party == "CD&V" & change == "fewer" & survey == "2023" ~ "Non-significant",
    party == "CD&V" & change == "fewer" & survey == "2024" ~ "Non-significant",
    
    party == "Open Vld" & change == "fewer" & survey == "2023" ~ "Non-significant",
    party == "Open Vld" & change == "fewer" & survey == "2024" ~ "Non-significant",
    
    party == "N-VA" & change == "fewer" & survey == "2023" ~ "Non-significant",
    party == "N-VA" & change == "fewer" & survey == "2024" ~ "Non-significant",
    
    party == "Vlaams Belang" & change == "fewer" & survey == "2023" ~ "Significant",
    party == "Vlaams Belang" & change == "fewer" & survey == "2024" ~ "Significant",
    
    party == "PVDA" & change == "same" & survey == "2023" ~ "Significant",
    party == "PVDA" & change == "same" & survey == "2024" ~ "Significant",
    
    party == "Groen" & change == "same" & survey == "2023" ~ "Non-significant",
    party == "Groen" & change == "same" & survey == "2024" ~ "Non-significant",
    
    party == "Vooruit" & change == "same" & survey == "2023" ~ "Significant",
    party == "Vooruit" & change == "same" & survey == "2024" ~ "Significant",
    
    party == "CD&V" & change == "same" & survey == "2023" ~ "Significant",
    party == "CD&V" & change == "same" & survey == "2024" ~ "Non-significant",
    
    party == "Open Vld" & change == "same" & survey == "2023" ~ "Non-significant",
    party == "Open Vld" & change == "same" & survey == "2024" ~ "Non-significant",
    
    party == "N-VA" & change == "same" & survey == "2023" ~ "Non-significant",
    party == "N-VA" & change == "same" & survey == "2024" ~ "Non-significant",
    
    party == "Vlaams Belang" & change == "same" & survey == "2023" ~ "Non-significant",
    party == "Vlaams Belang" & change == "same" & survey == "2024" ~ "Non-significant",
    
    party == "PVDA" & change == "more" & survey == "2023" ~ "Significant",
    party == "PVDA" & change == "more" & survey == "2024" ~ "Significant",
    
    party == "Groen" & change == "more" & survey == "2023" ~ "Non-significant",
    party == "Groen" & change == "more" & survey == "2024" ~ "Non-significant",
    
    party == "Vooruit" & change == "more" & survey == "2023" ~ "Non-significant",
    party == "Vooruit" & change == "more" & survey == "2024" ~ "Non-significant",
    
    party == "CD&V" & change == "more" & survey == "2023" ~ "Non-significant",
    party == "CD&V" & change == "more" & survey == "2024" ~ "Non-significant",
    
    party == "Open Vld" & change == "more" & survey == "2023" ~ "Non-significant",
    party == "Open Vld" & change == "more" & survey == "2024" ~ "Non-significant",
    
    party == "N-VA" & change == "more" & survey == "2023" ~ "Significant",
    party == "N-VA" & change == "more" & survey == "2024" ~ "Non-significant",
    
    party == "Vlaams Belang" & change == "more" & survey == "2023" ~ "Significant",
    party == "Vlaams Belang" & change == "more" & survey == "2024" ~ "Significant"))

# Assign numeric IDs for plotting order as above
prob.mass.df <- prob.mass.df %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  mutate(id = case_when(
    party == "PVDA" ~ 7,
    party == "Groen" ~ 6,
    party == "Vooruit" ~ 5,
    party == "CD&V" ~ 4,
    party == "Open Vld" ~ 3,
    party == "N-VA" ~ 2,
    party == "Vlaams Belang" ~ 1
  ))

# =============================================================================
# 48. Combine elite and mass data, cleanup, and prepare for plotting
# =============================================================================

# Combine elite and mass data into a single data frame for comparison
prob.diff.df <- rbind(prob.elite.df, prob.mass.df)

# Rename 'Vlaams Belang' party label to 'VB' for consistency/shortening
prob.diff.df <- prob.diff.df %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  mutate(party = recode(party, "Vlaams Belang" = "VB"))

# Reorder 'change' factor levels to ensure the desired facet order in plots
prob.diff.df$change <- factor(prob.diff.df$change, levels = c("fewer", "same", "more"))

# Create a formatted string label for the absolute difference (diff), rounded to one decimal place
prob.diff.df$diff_label <- sprintf("%.1f", prob.diff.df$diff)

# Pivot wider to separate 'P' and 'C' values
prob.diff.wide <- prob.diff.df %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  dplyr::select(party, change, survey, type, value, sig, diff_label) %>%
  pivot_wider(names_from = type, values_from = c(value, diff_label), names_sep = "_") %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  dplyr::mutate(
    show_diff_label = case_when(
      value_P > value_C ~ paste("Diff:", diff_label_P),
      value_C > value_P ~ paste("Diff:", diff_label_C),
      TRUE ~ NA_character_
    )
  )    

# Prepare long format for plotting
prob.diff.long <- prob.diff.wide %>%
  pivot_longer(cols = starts_with("value_"), names_to = "type", names_prefix = "value_", values_to = "value") %>%
  pivot_longer(cols = starts_with("diff_label_"), names_to = "diff_type", names_prefix = "diff_label_", values_to = "diff_label") %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  dplyr::filter(type == diff_type) %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  dplyr::select(-diff_type) %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  dplyr::mutate(id = row_number())

# Conditionally create a column for displaying difference labels ("Diff: x.x") on the plot
prob.diff.long <- prob.diff.long %>%
  dplyr::group_by(party, change, survey) %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  dplyr::mutate(
    # Show the diff label if the current group's value is greater than the opposite group (citizens vs politicians)
    show_diff_label = ifelse(
      (type == "P" & value > value[type == "C"]) | 
        (type == "C" & value > value[type == "P"]),
      paste("Diff:", diff_label),  # Include the "Diff:" prefix
      NA  # Otherwise leave blank (no label)
    )
  ) %>%
  ungroup()

prob.diff.long <- prob.diff.long %>%
  # --- Data transformation: ensure expected columns exist and types are correct (character vs factor vs numeric).
  dplyr::mutate(
    show_diff_label = ifelse(
      sig == "Significant" & !is.na(show_diff_label),
      paste0(show_diff_label, "*"),
      show_diff_label
    )
  )

prob.diff.long$party <- factor(prob.diff.long$party, 
                               levels = c("VB", "N-VA", "Open Vld", "CD&V", 
                                          "Vooruit", "Groen", "PVDA"))

# =============================================================================
# 49. Plot the difference in predicted probabilities
# =============================================================================

# Open a PNG device to save the plot with specified size and high resolution
png("prob_diff.png", units = "in", width = 8, height = 8, res = 1200)

# Create a ggplot object visualizing the gap in predicted probabilities across parties and change types
# --- Plotting: figure creation starts here. Review aesthetics (aes) and any facet/scales for readability.
ggplot(prob.diff.long, aes(x = party, y = value, group = party)) +
  geom_line() +  # Connect points within each party with lines
  geom_text(aes(label = type), color = "black", size = 3, show.legend = FALSE) +  # Add type label (C or P) near points
  geom_text(aes(label = show_diff_label), size = 3, vjust = -1.5, color = "black") +  # Conditionally display difference labels
  facet_grid(survey ~ change, scales = "free_y",  # Separate plots by change type with free y-axis scales
             labeller = labeller(change = c(fewer = "Fewer Seats", same = "Same as Now", more = "More Seats"))) +  # Custom facet titles
  theme_bw() +  # Apply a clean, white background theme
  labs(x = "Party", y = "Gap in Predicted Probabilities") +  # Axis labels
  scale_y_continuous(limits = c(0, 60), breaks = seq(0, 60, by = 10), expand = expansion(mult = c(0.15, 0.15))) +  # Y-axis limits and tick marks
  theme(
    axis.text.x = element_text(size = 10, hjust = .5),  # X-axis text formatting
    axis.text.y = element_text(size = 10),
    axis.title.x = element_text(margin = margin(t = 10)),  # Padding below x-axis title
    axis.title.y = element_blank(),  # Remove y-axis title
    strip.background = element_blank(),  # Remove facet strip background color
    strip.text = element_text(size = 11),  # Facet label font size
    panel.grid.minor = element_blank(),
    #legend.position = "bottom",  # Legend placement at the bottom
    legend.text = element_text(size = 10),
    legend.title = element_blank(),
    legend.position = "none",
    plot.title = element_blank()  # Remove any plot title
  ) +
  coord_flip()  # Flip axes to have parties on the y-axis and values on the x-axis

# Close the PNG device to finalize the saved image file
dev.off()

# =============================================================================
# 50. Forecast interaction: N-VA
# =============================================================================

# Seat forecast: N-VA

model.seats_nva <- clm(seats_nva ~ party_nva*type*survey + 
                         sex + age, family="binomial", data = flanders.data)

summary(model.seats_nva)

# Predicted probabilities of forecast for supporters and non-supporters among citizens and politicians

nva.fct.pred.fct.int.prop <- ggeffect(model.seats_nva, terms = c("type", "party_nva")) #categorical variables set to their proportion
nva.fct.pred.fct.int.prop <- ggeffect(model.seats_nva, terms = c("party_nva", "type")) #categorical variables set to their proportion

nva.fct.pred.fct.int.prop <- ggeffect(model.seats_nva, terms = c("type", "party_nva", "survey")) #categorical variables set to their proportion
nva.fct.pred.fct.int.prop <- ggeffect(model.seats_nva, terms = c("party_nva", "type", "survey")) #categorical variables set to their proportion


# Predicted probabilities of forecast for supporters and non-supporters

nva.fct.pred.party.prop <- ggeffect(model.seats_nva, terms = c("party_nva"))
nva.fct.pred.party.ref <- ggpredict(model.seats_nva, terms = c("party_nva"))

# Predicted probabilities of forecast for citizens and politicians

nva.fct.pred.type.prop <- ggeffect(model.seats_nva, terms = c("type"))
nva.fct.pred.type.ref <- ggpredict(model.seats_nva, terms = c("type"))

# =============================================================================
# 51. Forecast interaction: Groen
# =============================================================================

# Seat forecast: Groen

model.seats_groen <- clm(seats_groen ~ party_groen*type*survey + 
                           sex + age, family="binomial", data = flanders.data)

summary(model.seats_groen)

# Predicted probabilities of forecast for supporters and non-supporters among citizens and politicians

groen.fct.pred.fct.int.prop <- ggeffect(model.seats_groen, terms = c("type", "party_groen", "survey")) #categorical variables set to their proportion

groen.fct.pred.fct.int.prop <- ggeffect(model.seats_groen, terms = c("party_groen", "type")) #categorical variables set to their proportion

# Predicted probabilities of forecast for supporters and non-supporters

groen.fct.pred.party.prop <- ggeffect(model.seats_groen, terms = c("party_groen"))
groen.fct.pred.party.ref <- ggpredict(model.seats_groen, terms = c("party_groen"))

# Predicted probabilities of forecast for citizens and politicians

groen.fct.pred.type.prop <- ggeffect(model.seats_groen, terms = c("type"))
groen.fct.pred.type.ref <- ggpredict(model.seats_groen, terms = c("type"))

# =============================================================================
# 52. Forecast interaction: Vooruit
# =============================================================================

# Seat forecast: Vooruit

model.seats_vooruit <- clm(seats_vooruit ~ party_vooruit*type*survey + 
                             sex + age, family="binomial", data = flanders.data)

summary(model.seats_vooruit)

# Predicted probabilities of forecast for supporters and non-supporters among citizens and politicians

vooruit.fct.pred.fct.int.prop <- ggeffect(model.seats_vooruit, terms = c("type", "party_vooruit", "survey")) #categorical variables set to their proportion
vooruit.fct.pred.fct.int.prop <- ggeffect(model.seats_vooruit, terms = c("party_vooruit", "type")) #categorical variables set to their proportion

# Predicted probabilities of forecast for supporters and non-supporters

vooruit.fct.pred.party.prop <- ggeffect(model.seats_vooruit, terms = c("party_vooruit"))

# Predicted probabilities of forecast for citizens and politicians

vooruit.fct.pred.type.prop <- ggeffect(model.seats_vooruit, terms = c("type"))
vooruit.fct.pred.type.ref <- ggpredict(model.seats_vooruit, terms = c("type"))

# =============================================================================
# 53. Forecast interaction: CD&V
# =============================================================================

# Seat forecast: CD&V

model.seats_cdv <- clm(seats_cdv ~ party_cdv*type*survey + 
                         sex + age, family="binomial", data = flanders.data)

summary(model.seats_cdv)

# Predicted probabilities of forecast for supporters and non-supporters among citizens and politicians

cdv.fct.pred.fct.int.prop <- ggeffect(model.seats_cdv, terms = c("type", "party_cdv", "survey")) #categorical variables set to their proportion
cdv.fct.pred.fct.int.prop <- ggeffect(model.seats_cdv, terms = c("party_cdv", "type")) #categorical variables set to their proportion

# Predicted probabilities of forecast for supporters and non-supporters

cdv.fct.pred.party.prop <- ggeffect(model.seats_cdv, terms = c("party_cdv"))
cdv.fct.pred.party.ref <- ggpredict(model.seats_cdv, terms = c("party_cdv"))

# Predicted probabilities of forecast for citizens and politicians

cdv.fct.pred.type.prop <- ggeffect(model.seats_cdv, terms = c("type"))
cdv.fct.pred.type.ref <- ggpredict(model.seats_cdv, terms = c("type"))

# =============================================================================
# 54. Forecast interaction: Open Vld
# =============================================================================

# Seat forecast: Open Vld

model.seats_ovld <- clm(seats_ovld ~ party_ovld*type*survey + 
                          sex + age, family="binomial", data = flanders.data)

summary(model.seats_ovld)

# Predicted probabilities of forecast for supporters and non-supporters among citizens and politicians

ovld.fct.pred.fct.int.prop <- ggeffect(model.seats_ovld, terms = c("type", "party_ovld", "survey")) #categorical variables set to their proportion
ovld.fct.pred.fct.int.prop <- ggeffect(model.seats_ovld, terms = c("party_ovld", "type")) #categorical variables set to their proportion

# Predicted probabilities of forecast for supporters and non-supporters

ovld.fct.pred.party.prop <- ggeffect(model.seats_ovld, terms = c("party_ovld"))
ovld.fct.pred.party.ref <- ggpredict(model.seats_ovld, terms = c("party_ovld"))

# Predicted probabilities of forecast for citizens and politicians

ovld.fct.pred.type.prop <- ggeffect(model.seats_ovld, terms = c("type"))
ovld.fct.pred.type.ref <- ggpredict(model.seats_ovld, terms = c("type"))

# =============================================================================
# 55. Seat forecast: PVDA
# =============================================================================

# Seat forecast: PVDA

model.seats_pvda <- clm(seats_pvda ~ party_pvda*type*survey + 
                          sex + age, family="binomial", data = flanders.data)

summary(model.seats_pvda)

pvda.fct.pred.fct.int.prop <- ggeffect(model.seats_pvda, terms = c("party_pvda", "type", "survey")) #categorical variables set to their proportion

# Predicted probabilities of forecast for supporters and non-supporters

pvda.fct.pred.party.prop <- ggeffect(model.seats_pvda, terms = c("party_pvda"))
#pvda.fct.pred.party.ref <- ggpredict(model.seats_pvda, terms = c("party_pvda"))

# Predicted probabilities of forecast for citizens and politicians

pvda.fct.pred.type.prop <- ggeffect(model.seats_pvda, terms = c("type"))

# =============================================================================
# 56. Seat forecast: Vlaams Belang
# =============================================================================

# Seat forecast: Vlaams Belang

model.seats_vb <- clm(seats_vb ~ party_vb*type*survey + 
                        sex + age, family="binomial", data = flanders.data)

summary(model.seats_vb)

# Predicted probabilities of forecast for supporters and non-supporters

vb.fct.pred.party.prop <- ggeffect(model.seats_vb, terms = c("party_vb"))
#vb.fct.pred.party.ref <- ggpredict(model.seats_vb, terms = c("party_vb"))

# Predicted probabilities of forecast for citizens and politicians

vb.fct.pred.type.prop <- ggeffect(model.seats_vb, terms = c("type"))

# =============================================================================
# 57. Government forecast: PVDA
# =============================================================================

model.govt_pvda <- lm(govt_pvda ~ party_pvda*type*survey + sex + age, data = flanders.data)

# Linear prediction for supporters and non-supporters

pvda.fct.lpred.party.ref <- ggpredict(model.govt_pvda, terms = c("party_pvda"))

# Linear prediction for citizens and politicians

pvda.fct.lpred.type.ref <- ggpredict(model.govt_pvda, terms = c("type"))

# =============================================================================
# 58. Government forecast: Vooruit
# =============================================================================

model.govt_vooruit <- lm(govt_vooruit ~ party_vooruit*type*survey + sex + age, data = flanders.data)

# Linear prediction for supporters and non-supporters

vooruit.fct.lpred.party.ref <- ggpredict(model.govt_vooruit, terms = c("party_vooruit"))

# Linear prediction for citizens and politicians

vooruit.fct.lpred.type.ref <- ggpredict(model.govt_vooruit, terms = c("type"))

# =============================================================================
# 59. Government forecast: Groen
# =============================================================================

model.govt_groen <- lm(govt_groen ~ party_groen*type*survey + sex + age, data = flanders.data)

# Linear prediction for supporters and non-supporters

groen.fct.lpred.party.ref <- ggpredict(model.govt_groen, terms = c("party_groen"))

# Linear prediction for citizens and politicians

groen.fct.lpred.type.ref <- ggpredict(model.govt_groen, terms = c("type"))

# =============================================================================
# 60. Government forecast: CD&V
# =============================================================================

model.govt_cdv <- lm(govt_cdv ~ party_cdv*type*survey + sex + age, data = flanders.data)

# Linear prediction for supporters and non-supporters

cdv.fct.lpred.party.ref <- ggpredict(model.govt_cdv, terms = c("party_cdv"))

# Linear prediction for citizens and politicians

cdv.fct.lpred.type.ref <- ggpredict(model.govt_cdv, terms = c("type"))

# =============================================================================
# 61. Government forecast: Open Vld
# =============================================================================

model.govt_ovld <- lm(govt_ovld ~ party_ovld*type*survey + sex + age, data = flanders.data)

# Linear prediction for supporters and non-supporters

ovld.fct.lpred.party.ref <- ggpredict(model.govt_ovld, terms = c("party_ovld"))

# Linear prediction for citizens and politicians

ovld.fct.lpred.type.ref <- ggpredict(model.govt_ovld, terms = c("type"))

# =============================================================================
# 62. Government forecast: N-VA
# =============================================================================

model.govt_nva <- lm(govt_nva ~ party_nva*type*survey + sex + age, data = flanders.data)

# Linear prediction for supporters and non-supporters

nva.fct.lpred.party.ref <- ggpredict(model.govt_nva, terms = c("party_nva"))

# Linear prediction for citizens and politicians

nva.fct.lpred.type.ref <- ggpredict(model.govt_nva, terms = c("type"))

# =============================================================================
# 63. Government forecast: Vlaams Belang
# =============================================================================

model.govt_vb <- lm(govt_vb ~ party_vb*type*survey + sex + age, data = flanders.data)

# Linear prediction for supporters and non-supporters

vb.fct.lpred.party.ref <- ggpredict(model.govt_vb, terms = c("party_vb"))

# Linear prediction for citizens and politicians

vb.fct.lpred.type.ref <- ggpredict(model.govt_vb, terms = c("type"))

# =============================================================================
# 64. Government expectations: party-specific brier models + predicted values
# =============================================================================
# Each model predicts a party-specific Brier score outcome using:
#   party_support * type * survey  +  sex + age
#
# The party_support term varies by party (party_pvda, party_groen, ...).
# The 3-way interaction allows:
#   - the supporter vs non-supporter gap to vary by respondent type and year,
#   - the citizen vs politician gap to vary by supporter status and year,
#   - and all gaps to vary across survey waves.
# =============================================================================

# --- PVDA (government Brier) --------------------------------------------------
model.govt_pvda_brier <- lm(
  govt_pvda_brier ~ party_pvda * type * survey + sex + age,
  data = flanders.data
)

# Predicted values across supporter status only
pvda.brier.lpred.party.ref <- ggpredict(
  model.govt_pvda_brier,
  terms = c("party_pvda")
)

# Predicted values across respondent type only
pvda.brier.lpred.type.ref <- ggpredict(
  model.govt_pvda_brier,
  terms = c("type")
)

# --- Groen (government Brier) -------------------------------------------------
model.govt_groen_brier <- lm(
  govt_groen_brier ~ party_groen * type * survey + sex + age,
  data = flanders.data
)

groen.brier.lpred.party.ref <- ggpredict(
  model.govt_groen_brier,
  terms = c("party_groen")
)

groen.brier.lpred.type.ref <- ggpredict(
  model.govt_groen_brier,
  terms = c("type")
)

# --- Vooruit (government Brier) ----------------------------------------------
model.govt_vooruit_brier <- lm(
  govt_vooruit_brier ~ party_vooruit * type * survey + sex + age,
  data = flanders.data
)

vooruit.brier.lpred.party.ref <- ggpredict(
  model.govt_vooruit_brier,
  terms = c("party_vooruit")
)

vooruit.brier.lpred.type.ref <- ggpredict(
  model.govt_vooruit_brier,
  terms = c("type")
)

# --- CD&V (government Brier) -------------------------------------------------
model.govt_cdv_brier <- lm(
  govt_cdv_brier ~ party_cdv * type * survey + sex + age,
  data = flanders.data
)

cdv.brier.lpred.party.ref <- ggpredict(
  model.govt_cdv_brier,
  terms = c("party_cdv")
)

cdv.brier.lpred.type.ref <- ggpredict(
  model.govt_cdv_brier,
  terms = c("type")
)

# --- Open Vld (government Brier) ---------------------------------------------
model.govt_ovld_brier <- lm(
  govt_ovld_brier ~ party_ovld * type * survey + sex + age,
  data = flanders.data
)

ovld.brier.lpred.party.ref <- ggpredict(
  model.govt_ovld_brier,
  terms = c("party_ovld")
)

ovld.brier.lpred.type.ref <- ggpredict(
  model.govt_ovld_brier,
  terms = c("type")
)

# --- N-VA (government Brier) --------------------------------------------------
model.govt_nva_brier <- lm(
  govt_nva_brier ~ party_nva * type * survey + sex + age,
  data = flanders.data
)

nva.brier.lpred.party.ref <- ggpredict(
  model.govt_nva_brier,
  terms = c("party_nva")
)

nva.brier.lpred.type.ref <- ggpredict(
  model.govt_nva_brier,
  terms = c("type")
)

# --- Vlaams Belang (government Brier) ----------------------------------------
model.govt_vb_brier <- lm(
  govt_vb_brier ~ party_vb * type * survey + sex + age,
  data = flanders.data
)

vb.brier.lpred.party.ref <- ggpredict(
  model.govt_vb_brier,
  terms = c("party_vb")
)

vb.brier.lpred.type.ref <- ggpredict(
  model.govt_vb_brier,
  terms = c("type")
)

# =============================================================================
# 65. Government expectations: store model objects (parallel naming)
# =============================================================================
model.govt_pvda    <- lm(govt_pvda_brier    ~ party_pvda    * type * survey + sex + age, data = flanders.data)
model.govt_groen   <- lm(govt_groen_brier   ~ party_groen   * type * survey + sex + age, data = flanders.data)
model.govt_vooruit <- lm(govt_vooruit_brier ~ party_vooruit * type * survey + sex + age, data = flanders.data)
model.govt_cdv     <- lm(govt_cdv_brier     ~ party_cdv     * type * survey + sex + age, data = flanders.data)
model.govt_ovld    <- lm(govt_ovld_brier    ~ party_ovld    * type * survey + sex + age, data = flanders.data)
model.govt_nva     <- lm(govt_nva_brier     ~ party_nva     * type * survey + sex + age, data = flanders.data)
model.govt_vb      <- lm(govt_vb_brier      ~ party_vb      * type * survey + sex + age, data = flanders.data)

# =============================================================================
# 66. Seat expectations: contrast extraction across fit types (fewer/same/more)
# =============================================================================
# This section assumes you have, for each party, a precomputed data frame called:
#   model.seats.<party>.int
# that contains predicted values (e.g., from ggpredict or custom prediction code).
#
# Expected structure of model.seats.<party>.int:
#   - party_<party> : 0/1 supporter indicator column (e.g., party_pvda)
#   - type         : 0/1 respondent type (citizen vs politician)
#   - survey       : year (e.g., 2023, 2024)
#   - fit.More, fit.Same, fit.Fewer : predicted outcomes for each response option
#
# Goal:
#   For each party and each fit type, compute a set of contrasts capturing:
#     (a) supporter vs non-supporter differences within type/year
#     (b) politician vs citizen differences within supporter status/year
#     (c) changes across years for those contrasts
# =============================================================================

# --- Party ID map (stable ordering for tables/plots) --------------------------
party_list <- c(pvda = 1, groen = 2, vooruit = 3, cdv = 4, ovld = 5, nva = 6, vb = 7)

# --- Helper: retrieve a single predicted value from a model dataframe ---------
get_fit <- function(party_val, type_val, survey_val, fit_column, model_name) {
  # model_name is a string, e.g. "model.seats.pvda.int"
  # get() retrieves the object from the environment by name.
  model <- get(model_name)
  
  # IMPORTANT:
  # The filter below relies on a variable named `party` (loop variable) that
  # exists in the calling environment (the for-loop). This is deliberate in the
  # original code; we keep it and document it.
  #
  # It dynamically builds the correct supporter indicator column name, e.g.:
  #   "party_" + "pvda"  ->  "party_pvda"
  #
  # !!sym(...) converts the string column name to a symbol for dplyr to use.
  model %>%
    # --- Subset to the exact cell of interest (supporter status × type × year)
    filter(
      !!sym(paste0("party_", party)) == party_val,
      type   == type_val,
      survey == survey_val
    ) %>%
    # --- Pull the predicted value for the desired outcome column (fit type)
    pull(!!sym(fit_column))
}

# --- Fit “outcome columns” to process (three-category seat expectation) -------
fit_types <- c("fit.More", "fit.Same", "fit.Fewer")

# Container for results across fit_types
all_results <- list()

# Loop over each outcome (“More seats”, “Same”, “Fewer seats”)
for (fit_type in fit_types) {
  fit_results <- list()
  
  # Loop over each party and compute the full set of contrasts for that party
  for (party in names(party_list)) {
    # Expected object name: model.seats.pvda.int, model.seats.groen.int, ...
    model_name <- paste0("model.seats.", party, ".int")
    
    # -------------------------------------------------------------------------
    # Contrast definitions (all computed on the predicted-value scale)
    #
    # Notation:
    #   party_val: 1=supporter, 0=non-supporter
    #   type_val : 0=citizens, 1=politicians   (assumed from the code)
    #   survey   : 2023 vs 2024
    #
    # Each contrast is defined as a difference of two predicted values.
    # Multiplying by 100 later expresses as percentage-point differences.
    # -------------------------------------------------------------------------
    diffs <- list(
      # Supporters vs Non-supporters among Citizens (type=0), within year
      cf_pid_2023 = get_fit(1, 0, 2023, fit_type, model_name) - get_fit(0, 0, 2023, fit_type, model_name),
      cf_pid_2024 = get_fit(1, 0, 2024, fit_type, model_name) - get_fit(0, 0, 2024, fit_type, model_name),
      
      # Supporters vs Non-supporters among Politicians (type=1), within year
      pf_pid_2023 = get_fit(1, 1, 2023, fit_type, model_name) - get_fit(0, 1, 2023, fit_type, model_name),
      pf_pid_2024 = get_fit(1, 1, 2024, fit_type, model_name) - get_fit(0, 1, 2024, fit_type, model_name),
      
      # Politicians vs Citizens among Non-supporters (party=0), within year
      type_ns_2023 = get_fit(0, 1, 2023, fit_type, model_name) - get_fit(0, 0, 2023, fit_type, model_name),
      type_ns_2024 = get_fit(0, 1, 2024, fit_type, model_name) - get_fit(0, 0, 2024, fit_type, model_name),
      
      # Politicians vs Citizens among Supporters (party=1), within year
      type_s_2023 = get_fit(1, 1, 2023, fit_type, model_name) - get_fit(1, 0, 2023, fit_type, model_name),
      type_s_2024 = get_fit(1, 1, 2024, fit_type, model_name) - get_fit(1, 0, 2024, fit_type, model_name),
      
      # Across-year differences for Citizens (supporters in 2024 vs non-supporters in 2023)
      # NOTE: this is not a within-group year change; it mixes supporter statuses across years,
      # exactly as written in your original code.
      cf_year_pid = get_fit(1, 0, 2024, fit_type, model_name) - get_fit(0, 0, 2023, fit_type, model_name),
      
      # Across-year differences for Politicians (supporters in 2024 vs non-supporters in 2023)
      # Same note as above re: mixed statuses.
      pf_year_pid = get_fit(1, 1, 2024, fit_type, model_name) - get_fit(0, 1, 2023, fit_type, model_name),
      
      # Across-year differences for Non-supporters (politicians 2024 vs citizens 2023)
      year_type_ns = get_fit(0, 1, 2024, fit_type, model_name) - get_fit(0, 0, 2023, fit_type, model_name),
      
      # Across-year differences for Supporters (politicians 2024 vs citizens 2023)
      year_type_s  = get_fit(1, 1, 2024, fit_type, model_name) - get_fit(1, 0, 2023, fit_type, model_name)
    )
    
    # Store results in long format for easy faceting/plotting/export
    fit_results[[party]] <- data.frame(
      id       = party_list[party],
      party    = party,
      metric   = names(diffs),
      value    = round(unlist(diffs) * 100, 2),  # convert to percentage points
      fit_type = fit_type,
      row.names = NULL
    )
  }
  
  # Bind parties for this fit_type
  all_results[[fit_type]] <- bind_rows(fit_results)
}

# --- Final tidy data for seats contrasts -------------------------------------
final_df <- bind_rows(all_results)

# Order and label fit_type for plotting (Fewer / Same / More)
final_df$fit_type <- factor(
  final_df$fit_type,
  levels = c("fit.Fewer", "fit.Same", "fit.More"),
  labels = c("Fewer", "Same", "More")
)

# Order and label parties for plotting/legend
final_df$party <- factor(
  final_df$party,
  levels = c("pvda", "groen", "vooruit", "cdv", "ovld", "nva", "vb"),
  labels = c("PVDA", "Groen", "Vooruit", "CD&V", "Open Vld", "N-VA", "Vlaams Belang")
)

# Color palette for party fills (used in the plot)
party_colors <- c(
  "PVDA"          = "#990000",
  "Groen"         = "#04847A",
  "Vooruit"       = "#FF0000",
  "CD&V"          = "#FF6000",
  "Open Vld"      = "#005DAA",
  "N-VA"          = "#FCBD1B",
  "Vlaams Belang" = "black"
)

# Human-readable labels for each contrast “metric”
final_df$metric <- factor(
  final_df$metric,
  levels = c(
    "cf_pid_2023", "cf_pid_2024",
    "pf_pid_2023", "pf_pid_2024",
    "type_ns_2023", "type_ns_2024",
    "type_s_2023", "type_s_2024",
    "cf_year_pid", "pf_year_pid",
    "year_type_ns", "year_type_s"
  ),
  labels = c(
    "Citizens: Supporters v.\nNon Supporters, 2023",
    "Citizens: Supporters v.\nNon Supporters, 2024",
    "Politicians: Supporters v.\nNon Supporters, 2023",
    "Politicians: Supporters v.\nNon Supporters, 2024",
    "Non Supporters: Politicians v.\nCitizens, 2023",
    "Non Supporters: Politicians v.\nCitizens, 2024",
    "Supporters: Politicians v.\nCitizens, 2023",
    "Supporters: Politicians v.\nCitizens, 2024",
    "Citizens: Supporters v.\nNon Supporters, 2023–2024",
    "Politicians: Supporters v.\nNon Supporters, 2023–2024",
    "Non Supporters: Politicians v.\nCitizens, 2023–2024",
    "Supporters: Politicians v.\nCitizens, 2023–2024"
  )
)

# --- Export seats contrast table to Excel (used in paper, section 4.2) --------
write_xlsx(
  final_df,
  "party_differences_seats_all_fits.xlsx"
) # Results reported in section 4.2 of the paper (seat expectations)

# =============================================================================
# 67. Seat expectations: faceted bar plot (png export)
# =============================================================================
# Plot shows, for each contrast metric (facet), the difference (pp) by:
#   x = seat expectation outcome (Fewer/Same/More),
#   fill = party,
#   grouped via dodge (side-by-side party bars).
#
# Two vertical lines visually separate the three outcome columns.
# =============================================================================

# Positions between the three fit_type categories (Fewer | Same | More)
vlines_df <- data.frame(xintercept = c(1.5, 2.5))

# High-resolution PNG export for publication-quality rendering
png("diff_all_g.png", units = "in", width = 12, height = 10, res = 1200)

ggplot(final_df, aes(x = fit_type, y = value, fill = party)) +
  # --- Vertical separators between outcome categories ------------------------
geom_vline(
  data = vlines_df,
  aes(xintercept = xintercept),
  inherit.aes = FALSE,
  color = "gray50",
  linewidth = 0.4
) +
  # --- Bars: one bar per party per outcome category --------------------------
geom_bar(stat = "identity", position = "dodge") +
  # --- Small multiples by contrast metric ------------------------------------
facet_wrap(~ metric, scales = "free_y", ncol = 4) +
  # --- Party color mapping ----------------------------------------------------
scale_fill_manual(values = party_colors) +
  # --- Add small padding to y-range so labels/bars don’t collide --------------
scale_y_continuous(expand = expansion(mult = c(0.05, 0.05))) +
  # --- Axis/legend labels -----------------------------------------------------
labs(
  x = "More or Less Seats?",
  y = "Difference (percentage points)",
  fill = "Party"
) +
  theme_minimal() +
  theme(
    panel.border    = element_rect(color = "black", fill = NA, linewidth = 0.6),
    strip.text      = element_text(size = 11),
    axis.text.x     = element_text(size = 11),
    axis.text.y     = element_text(size = 11),
    axis.title.x    = element_text(margin = margin(t = 12)),
    axis.title.y    = element_text(margin = margin(r = 15)),
    legend.position = "bottom",
    legend.margin   = margin(t = 12),
    plot.title      = element_blank()
  )

dev.off()

# =============================================================================
# 68. Government expectations: contrast table (single "fit" column per party)
# =============================================================================
# Expected structure of model.govt.<party>.df:
#   - party  : 0/1 supporter indicator (generic column name)
#   - type   : 0/1 respondent type
#   - survey : year (e.g., 2023, 2024)
#   - fit    : predicted value (single outcome column)
#
# Goal:
#   Compute the same suite of contrasts as above, but for the government
#   expectations outcome.
# =============================================================================

# --- Party ID map (duplicated in original; kept for clarity/standalone runs) --
party_list <- c(pvda = 1, groen = 2, vooruit = 3, cdv = 4, ovld = 5, nva = 6, vb = 7)

# --- Helper: retrieve a single predicted “fit” from a govt model dataframe ----
get_fit <- function(party_val, type_val, survey_val, model_name) {
  model <- get(model_name)
  
  model %>%
    # --- Subset to the exact cell of interest --------------------------------
  filter(
    party  == party_val,
    type   == type_val,
    survey == survey_val
  ) %>%
    # --- Pull the predicted value --------------------------------------------
  pull(fit)
}

all_results <- list()

# Loop over parties, compute contrasts, bind into one long tidy table
for (party in names(party_list)) {
  model_name <- paste0("model.govt.", party, ".df")
  
  diffs <- list(
    cf_pid_2023  = get_fit(1, 0, 2023, model_name) - get_fit(0, 0, 2023, model_name),
    cf_pid_2024  = get_fit(1, 0, 2024, model_name) - get_fit(0, 0, 2024, model_name),
    pf_pid_2023  = get_fit(1, 1, 2023, model_name) - get_fit(0, 1, 2023, model_name),
    pf_pid_2024  = get_fit(1, 1, 2024, model_name) - get_fit(0, 1, 2024, model_name),
    type_ns_2023 = get_fit(0, 1, 2023, model_name) - get_fit(0, 0, 2023, model_name),
    type_ns_2024 = get_fit(0, 1, 2024, model_name) - get_fit(0, 0, 2024, model_name),
    type_s_2023  = get_fit(1, 1, 2023, model_name) - get_fit(1, 0, 2023, model_name),
    type_s_2024  = get_fit(1, 1, 2024, model_name) - get_fit(1, 0, 2024, model_name),
    cf_year_pid  = get_fit(1, 0, 2024, model_name) - get_fit(0, 0, 2023, model_name),
    pf_year_pid  = get_fit(1, 1, 2024, model_name) - get_fit(0, 1, 2023, model_name),
    year_type_ns = get_fit(0, 1, 2024, model_name) - get_fit(0, 0, 2023, model_name),
    year_type_s  = get_fit(1, 1, 2024, model_name) - get_fit(1, 0, 2023, model_name)
  )
  
  all_results[[party]] <- data.frame(
    id     = party_list[party],
    party  = party,
    metric = names(diffs),
    value  = round(unlist(diffs), 2),  # NOTE: not multiplied by 100 in original
    row.names = NULL
  )
}

# Tidy long table across parties
final_govt_df <- bind_rows(all_results)

# Party labels for readability
final_govt_df$party <- factor(
  final_govt_df$party,
  levels = c("pvda", "groen", "vooruit", "cdv", "ovld", "nva", "vb"),
  labels = c("PVDA", "Groen", "Vooruit", "CD&V", "Open Vld", "N-VA", "Vlaams Belang")
)

# Metric labels (same as seats section)
final_govt_df$metric <- factor(
  final_govt_df$metric,
  levels = c(
    "cf_pid_2023", "cf_pid_2024",
    "pf_pid_2023", "pf_pid_2024",
    "type_ns_2023", "type_ns_2024",
    "type_s_2023", "type_s_2024",
    "cf_year_pid", "pf_year_pid",
    "year_type_ns", "year_type_s"
  ),
  labels = c(
    "Citizens: Supporters v.\nNon Supporters, 2023",
    "Citizens: Supporters v.\nNon Supporters, 2024",
    "Politicians: Supporters v.\nNon Supporters, 2023",
    "Politicians: Supporters v.\nNon Supporters, 2024",
    "Non Supporters: Politicians v.\nCitizens, 2023",
    "Non Supporters: Politicians v.\nCitizens, 2024",
    "Supporters: Politicians v.\nCitizens, 2023",
    "Supporters: Politicians v.\nCitizens, 2024",
    "Citizens: Supporters v.\nNon Supporters, 2023–2024",
    "Politicians: Supporters v.\nNon Supporters, 2023–2024",
    "Non Supporters: Politicians v.\nCitizens, 2023–2024",
    "Supporters: Politicians v.\nCitizens, 2023–2024"
  )
)

# --- Export government contrast table to Excel (paper section 4.2) ------------
write_xlsx(
  final_govt_df,
  "party_differences_govt_all_fits.xlsx"
) # Results reported in section 4.2 of the paper (government expectations)

# =============================================================================
# END OF SCRIPT
# =============================================================================